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Abstract 

Sound  propagation  in  shallow  water  is  highly  dependent  on  the  interaction  of  the 
sound  field  with  the  bottom.  In  order  to  fully  understand  this  problem,  it  is  neces¬ 
sary  to  obtain  reliable  estimates  of  bottom  geoacoustic  properties  that  can  be  used 
in  acoustic  propagation  codes.  In  this  thesis,  perturbative  inversion  methods  and 
exact  inverse  methods  are  discussed  as  a  means  for  inferring  geoacoustic  properties 
of  the  bottom.  For  each  of  these  methods,  the  input  data  to  the  inversion  is  the 
horizontal  wavenumber  spectrum  of  a  point-source  acoustic  field.  The  main  thrust 
of  the  thesis  work  concerns  extracting  horizontal  wavenumber  content  for  fully  three- 
dimensionally  varying  waveguide  environments.  In  this  context,  a  high-resolution 
autoregressive  (AR)  spectral  estimator  was  applied  to  determine  wavenumber  con¬ 
tent  for  short  aperture  data.  As  part  of  this  work,  the  AR  estimator  was  examined 
for  its  ability  to  detect  discrete  wavenumbers  in  the  presence  of  noise  and  also  to 
resolve  closely  spaced  wavenumbers  for  short  aperture  data.  As  part  of  a  geoacous¬ 
tic  inversion  workshop,  the  estimator  was  applied  to  extract  horizontal  wavenumber 
content  for  synthetic  pressure  field  data  with  range- varying  geoacoustic  properties  in 
the  sediment.  The  resulting  wavenumber  content  was  used  as  input  data  to  a  per¬ 
turbative  inverse  algorithm  to  determine  the  sound  speed  profile  in  the  sediment.  It 
was  shown  using  the  high-resolution  wavenumber  estimator  that  both  the  shape  and 
location  of  the  range- variability  in  the  sediment  could  be  determined.  The  estima¬ 
tor  was  also  applied  to  determine  wavenumbers  for  synthetic  data  where  the  water 
column  sound  speed  contained  temporal  variations  due  to  the  presence  of  internal 
waves.  It  was  shown  that  reliable  estimates  of  horizontal  wavenumbers  could  be  ob¬ 
tained  that  are  consistent  with  the  boundary  conditions  of  the  waveguide.  The  Modal 
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Mapping  Experiment  (MOM AX),  an  experimental  method  for  measuring  the  full  spa¬ 
tial  variability  of  a  propagating  sound  field  and  its  corresponding  modal  content  in 
two-dimensions,  is  also  discussed.  The  AR  estimator  is  applied  to  extract  modal  con¬ 
tent  from  the  real  data  and  interpreted  with  respect  to  source/receiver  motion  and 
geometry.  For  a  moving  source,  it  is  shown  that  the  wavenumber  content  is  Doppler 
shifted.  A  method  is  then  described  that  allows  the  direct  measure  of  modal  group 
velocities  from  Doppler  shifted  wavenumber  spectra.  Finally,  numerical  studies  are 
presented  addressing  the  practical  issues  associated  with  using  MOMAX  type  data 
in  the  exact  inversion  method  of  Gelfand-Levitan. 

Thesis  Supervisor:  George  V.  Frisk 
Title:  Senior  Scientist,  WHOI 
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Wavenumber  and  group  velocities  estimated  from  moving  source  spec¬ 
tra  compared  to  true  values.  Source  speed  was  10  m/s . 


Chapter  1 


Introduction  and  Overview 


1.1  Motivation 

Any  review  of  the  underwater  acoustic  or  ocean  science  literature  will  immediately 
impress  upon  the  reader  the  importance  of  acoustics  as  a  tool  for  understanding  not 
only  acoustic  problems  such  as  propagation  and  scattering,  but  many  other  ocean 
science  problems  as  well,  including  measurements  of  ocean  temperature  and  currents, 
monitoring  of  marine  mammal  activity,  and  characterization  of  different  fish  species, 
to  name  a  few.  Whatever  the  application,  each  benefits  from  increased  knowledge 
of  the  propagation  environment.  Consistent  with  this,  in  a  much  referenced  paper 
by  Hamilton[36],  a  case  was  made  for  the  importance  of  geoacoustic  modeling  of  the 
seabed.  Hamilton  asserts  that  knowledge  of  geoacoustic  models  is  basic  to  the  under¬ 
standing  of  underwater  acoustics.  His  statement  of  twenty  years  ago  is  particularly 
true  now  with  today’s  requirement  of  rapid  assessment  where  complex  models  based 
on  realistic  input  parameters  are  used  for  prediction  purposes.  In  shallow  water  envi¬ 
ronments,  where  the  sound  field  is  constantly  interacting  with  the  seafloor,  knowledge 
of  bottom  geoacoustic  properties  is  critical  to  making  accurate  predictions.  It  is  this 
scenario  that  calls  for  a  robust  means  of  determining  the  geoacoustic  properties  of 
the  seafloor. 
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Stratified  Waveguide  Model 


Figure  1-1:  Simple  environmental  model  representing  horizontally  stratified  waveg¬ 
uide 

1.2  Objectives 

The  objective  of  this  thesis  work  is  to  investigate  methods  for  inferring  the  properties 
of  the  seafloor  from  spectral  measurements  of  point-source  acoustic  pressure  field  data 
in  shallow-water  environments.  The  data  consist  of  complex  pressure  measured  on 
synthetic  aperture  horizontal  planar  arrays  for  a  point  source  at  constant  depth  emit¬ 
ting  several  discrete  tones  [31].  A  canonical  model  is  represented  in  Figure  1-1.  In 
this  context,  the  spectral  representation  of  the  field  is  given  by  the  depth-dependent 
Green’s  function  versus  horizontal  wavenumber.  As  will  be  shown,  the  basis  for  the 
geoacoustic  inverse  problem  is  the  dependence  of  the  Green’s  function  on  the  local 
properties  of  the  waveguide.  Theory  and  experiment  have  confirmed  this  relationship 
through  application  of  Hankel  transform-based  methods  for  extracting  the  wavenum- 
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ber  spectrum  for  data  measured  on  synthetic  aperture  line  arrays  [26]  [55].  The  work 
described  in  this  thesis  differs  from  previous  work  in  that  synthetic  apertures  are 
formed  along  arbitrary  horizontal  paths  between  source  and  receiver  positions  as  op¬ 
posed  to  being  formed  along  radials  at  constant  bearing.  These  measurements  are 
unique  in  that  they  use  precision  Global  Positioning  System  (GPS)  measurements  of 
source  and  receiver  buoy  positions  to  create  the  synthetic  aperture  used  for  processing. 
Examination  of  the  point  source  spectra  determined  for  sub-apertures  of  these  hori¬ 
zontal  arrays  yields  information  about  spatial  variability  in  the  seafloor.  Emphasis  is 
placed  on  extracting  modal  content  and  interpreting  modal  evolution  in  the  context 
of  the  measurement  geometry.  Where  appropriate,  this  includes  Doppler  effects  due 
to  source/receiver  motion  on  the  measured  spectra.  In  this  thesis,  a  procedure  for 
the  direct  measurement  of  modal  group  velocity  is  presented  through  measurement 
of  an  observed  Doppler  shift.  Using  this  result,  the  effect  of  Doppler  can  be  removed 
from  the  wavenumber  estimates.  Along  these  same  lines,  another  objective  of  this 
work  is  to  achieve  maximum  resolution  of  discrete  wavenumbers  given  short  aperture 
field  data  to  obtain  high  resolution  geospatial  maps  of  wavenumber  content.  In  this 
effort,  the  use  of  high  resolution  spectral  estimation  techniques  are  examined  as  well 
as  the  effect  of  sound  speed  variations  in  the  water  column  on  extraction  of  closely 
spaced  spectral  values. 

Given  estimates  of  the  point  source  spectrum,  the  next  objective  is  to  use  the 
spectral  data  as  input  for  the  geoacoustic  inverse  problem.  The  experiments  for 
this  work  were  designed  primarily  as  a  means  to  provide  input  data  to  an  inverse 
method  based  on  making  linear  perturbations  to  a  background  starting  model.  This 
approach  uses  the  discrete  portion  of  the  point  source  spectrum  given  by  peaks  in 
the  Green’s  functions  as  input  data  to  the  inverse  problem.  Emphasis  is  placed 
on  examining  the  applicability  of  the  measurement  data  obtained  on  the  arbitrary 
horizontal  arrays  to  the  perturbative  inverse  approach.  During  the  course  of  this  work, 
it  was  suggested  that  point-source  spectral  measurements  as  described  might  also  be 
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used  as  input  data  to  one  of  the  exact  methods  for  geoacoustic  inversion  that  does 
not  require  perturbations  or  iterations.  Thus,  a  second  inversion  approach  studied 
is  an  adaptation  of  the  Gelfand-Levitan  method  for  geoacoustic  inversion  using  the 
point-source  spectral  data  for  shallow  water.  In  this  approach  both  the  continuous 
and  discrete  portion  of  the  Green’s  function  are  used  to  obtain  the  input  data  for  the 
inverse  algorithm.  It  will  be  shown  that  from  these  measurements,  it  may  be  possible 
to  estimate  the  reflection  coefficient  at  the  water-bottom  interface.  The  input  for  the 
inversion  scheme  is  then  the  Fourier  transform  of  the  plane-wave  reflection  coefficient. 
Numerical  results  will  be  presented  for  this  case  with  recommendations  made  for 
future  experiments. 

Clearly,  by  combining  the  spectral  estimation  efforts  and  the  inversion  efforts,  the 
local  nature  of  the  Green’s  function  can  be  exploited  to  obtain  fully  three-dimensional 
geoacoustic  inversions. 

1.3  B  ackgr  ound 

The  history  of  seafloor  characterization  in  the  underwater  acoustics  community  is  a 
long  one.  Many  different  measurement  and  analysis  techniques  have  been  developed 
through  the  years  to  estimate  bathymetry,  seafloor  roughness,  and  sediment  compo¬ 
sition.  In  probing  the  seafloor  for  its  geoacoustic  properties,  several  major  classes  of 
approaches  have  emerged.  Frisk  [27]  breaks  these  classes  into  four  categories  that  in¬ 
clude  inversion  for  bottom  properties  based  on  specific  features,  iteration  of  forward 
models,  perturbative  inversion,  or  exact  inverse  methods.  Included  in  the  class  of 
iterating  forward  models  are  optimization  methods  such  as  simulated  annealing  or 
genetic  algorithms.  For  the  measurements  and  data  to  be  discussed  in  regard  to  this 
thesis  work,  any  of  the  latter  three  inversion  method  classes  could  be  appropriate. 
However,  of  particular  interest  are  the  perturbative  inverse  approach  and  exact  meth¬ 
ods.  It  should  be  noted  that  there  is  a  bit  of  overlap  in  these  classifications  in  that 
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both  the  perturbative  inverse  method  and  the  iteration  of  forward  models  involve 
iteratively  making  perturbations  to  a  starting  model  based  on  a  priori  information 
about  the  environment  combined  with  measured  field  data.  The  differences  arise  from 
the  nature  of  the  perturbations,  as  well  as  the  cost  function  defined  to  evaluate  the 
results.  These  differences  lead  to  a  substantial  difference  in  the  number  of  iterations 
required  and  are  discussed  further  in  the  following  section  where  both  approaches  are 
categorically  referred  to  as  iterative  methods.  For  any  of  these  iterative  methods, 
results  from  an  inversion  are  compared  to  measured  data  and  the  model  is  perturbed 
again  using  the  comparison  information.  Iteration  continues  until  the  model/data 
comparison  falls  within  an  acceptable  range  of  the  observed  values.  Exact  methods 
differ  from  optimization  methods  in  that  they  do  not  require  a  background  model 
for  determination  of  the  model  parameters.  These  methods  are  exact  in  the  sense 
that  they  have  either  analytic  or  numerical  forms  of  solution  which  make  direct  use 
of  the  field  data  in  the  reconstruction  of  the  geoacoustic  properties.  Consequently, 
these  methods  do  not  suffer  from  any  reliance  on  the  approximations  made  for  the 
forward-problem  model,  or  the  sometimes  non-unique  nature  of  the  perturbed  back¬ 
ground  model  solutions.  Following  is  a  brief  summary  of  previous  work  done  using 
the  two  types  of  methods. 

1.3.1  Iterative  Methods 

Much  of  the  recent  geoacoustic  inversion  literature  has  focused  on  algorithms  which 
seek  to  determine  a  set  of  model  parameters  which  minimize  the  mismatch  between 
measured  and  simulated  field  data.  This  problem  becomes  increasingly  difficult  and 
computationally  intensive  as  the  number  of  model  parameters  to  be  determined  is 
increased.  Further  complicating  the  problem  is  the  presence  of  local  minima  inher¬ 
ent  in  the  inverse  problem.  Rajan  et  al.  [62]  [66]  employed  several  different  iterative 
approaches  to  perturb  a  starting  model  using  measurements  of  the  complex  pres¬ 
sure  field  as  input  data.  Using  short  aperture  Hankel  transforms  to  obtain  eigenvalue 
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information,  Frisk  et  al.  [26]  used  a  perturbative-inverse  algorithm  to  obtain  compres- 
sional  wave  speed  profiles  for  two  adjacent,  but  geo-acoustically  distinct,  regions  in 
Nantucket  Sound.  Ohta  and  Frisk[54][55]  used  a  short-window  sliding  transform  and 
iteration  of  forward  models  to  further  extend  these  methods  to  the  range-dependent 
case.  The  perturbative  inversion  methods  used  by  these  authors  were  based  on  linear 
or  quasi-linear  inverse  algorithms.  To  ensure  linearity,  the  parameter  search  space 
was  limited  to  small  perturbations  around  the  background  model.  Further,  the  op¬ 
erator  on  the  parameter  that  gives  an  estimate  of  the  data  is  assumed  to  be  either 
monotonic,  or  Frechet  differentiable,  a  condition  that  allows  the  linearization  of  a 
nonlinear  function  [64].  These  requirements  ensure  that  the  inversion  algorithm  trav¬ 
els  downhill  and  will  arrive  at  the  global  minimum  for  the  cost  function  as  determined 
by  linearizing  the  problem  around  the  assumed  background  model.  The  benefits  of 
the  perturbative  inverse  method  include  the  following:  the  search  parameter  space 
derived  from  the  data  limits  the  results  to  realistic  values  for  the  seafloor;  the  linear 
inverse  equations  allows  for  the  determination  of  the  variance  and  resolution  of  the 
parameter  estimates;  and  the  results  do  not  depend  on  and  do  not  require  multiple 
full  forward  modeling  runs.  A  limitation  of  these  methods  are  their  sensitivity  to 
the  starting  model  and  the  potential  for  arriving  at  at  a  non-unique  solution  due 
to  minimizing  the  cost  function  based  on  an  incorrect  starting  background  model. 
To  overcome  these  shortcomings,  other  nonlinear  approaches  have  been  investigated 
which  allow  a  more  thorough  search  of  the  model  parameter  space  by  employing  ran¬ 
dom  processes  to  perturb  the  starting  model.  The  random  nature  of  the  perturbation 
makes  them  less  sensitive  to  the  starting  model  and  avoids  trapping  in  the  local  min¬ 
ima  along  with  removing  the  monotonic  requirement  on  the  data.  However,  because 
of  their  random  nature,  global  methods  tend  to  converge  more  slowly  and  are  more 
computationally  expensive.  Further,  because  global  methods  require  repeated  appli¬ 
cation  of  forward  modeling  runs,  they  are  dependent  on  the  forward  model  used.  Two 
of  the  more  popular  global  methods  that  have  been  used  successfully  in  geoacoustic 
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inversions  are  simulated  annealing  (SA)  and  genetic  algorithms  (GA).  SA  is  based  on 
the  thermodynamic  annealing  principles  governing  crystal  growth  and  seeks  a  mini¬ 
mum  in  parameter  space  through  random  perturbations  of  the  model  with  acceptance 
of  uphill  steps  probabilistically  decreasing  as  time  increases.  Kuperman  et  al.  [45] 
devised  a  fast  SA  algorithm  to  invert  for  wave-speeds  and  attenuation  for  slightly 
range-dependent  environments.  Collins  et  al.  [15]  used  simulated  annealing  based  on 
a  full-field  parabolic  equation  model  to  invert  for  range-  and  depth-dependent  bottom 
geoacoustic  parameters.  GA  are  analogous  to  genetic  evolution  where  perturbations 
to  a  model  population  are  related  to  genetic  crossover,  or  recombination,  and  muta¬ 
tions  resulting  in  a  new  model  of  a  given  fitness.  Models  with  high  fitness,  indicating 
a  low  mismatch  between  parent  and  offspring,  are  kept,  while  low  fitness  models  are 
rejected.  Different  probabilistic  models  can  be  introduced  into  the  algorithm  to  fur¬ 
ther  refine  the  fitness  criterion  for  the  selection  of  parents  in  the  formation  of  new 
offspring.  GA  was  first  used  on  synthetic  data  by  Gerstoft[35],  and  then  used  on  real 
broadband  acoustic  data  by  Hermand  and  Gerstoft  [39]  to  determine  bottom  sedi¬ 
ment  properties  and  layer  thicknesses.  Combining  the  quick  convergence  properties 
of  the  local  approaches  with  the  full  parameter  space  searches  of  global  methods, 
Fallat  and  Dosso[23]  have  devised  a  hybrid  approach  based  on  downhill  simplex  and 
simulated  annealing  to  invert  for  geoacoustic  properties.  Also  seeking  to  minimize 
the  mismatch  between  measured  and  simulated  data,  their  method  was  shown  to  con¬ 
verge  demonstrably  faster  than  either  a  local  or  global  method  alone  for  a  test  case 
where  the  geoacoustic  parameters  were  known. 

Each  of  the  inversion  methods  described  above  has  its  advantages  and  disadvan¬ 
tages.  The  computational  expense  of  the  various  methods  will  be  addressed  in  a 
forthcoming  report  from  a  recent  workshop  on  geoacoustic  inversion  techniques.  The 
report  should  provide  a  good  summary  of  the  computational  requirements  and  effi¬ 
ciency  of  various  inverse  methods  [13].  Putting  the  other  differences  aside,  the  work 
is  this  thesis  was  motivated  by  experimental  work  designed  to  provide  input  data  for 


23 


a  perturbative  inverse  method.  Therefore,  much  of  the  development  in  this  thesis  is 
placed  in  the  context  of  perturbative  inversion. 

1.3.2  Exact  Inverse  Methods 

Much  of  the  exact  inverse  work  in  determining  geoacoustic  properties  considered  for 
this  thesis  is  found  in  the  geophysical  literature  and  is  summarized  by  Merab[50].  The 
geophysical  problem  of  interest  is  largely  based  on  physics  and  mathematics  in  the 
form  of  the  solution  to  inverse  problems  of  quantum  scattering  theory.  The  connection 
between  scattering  theory  and  the  acoustic  problem  is  made  by  showing  the  geoacous¬ 
tic  inverse  problem  to  be  equivalent  to  determining  the  potential  in  one-dimensional 
inverse  scattering  problems  for  the  Schrodinger  equation.  Of  particular  relevance  to 
the  solution  of  this  problem  are  the  methods  of  Gelfand  and  Levitan[34],  and  the 
trace  methods  formulated  by  Deift  and  Trubowitz[17]  (DT).  The  formulation  of  these 
two  solution  alternatives  begin  with  equivalent  equations,  boundary  conditions,  and 
assumptions. 

The  major  difference  between  the  two  methods  of  solution  is  in  the  choice  of  input 
data  for  the  inversion.  The  DT  trace  method  uses  the  plane- wave  reflection  coefficient 
as  a  function  of  wavenumber  as  input  data  for  the  simultaneous  solution  of  a  coupled 
pair  of  equations.  The  solution  is  given  in  terms  of  a  potential  q{z)  that  is  directly 
related  to  the  reflection  coefficient.  It  will  be  shown  in  chapter  2  that  the  plane-wave 
reflection  coefficient  is  directly  related  to  the  sediment  properties  thus  establishing 
the  relationship  between  the  potential  and  the  sound  speed  profile  in  the  bottom.  In 
its  most  general  form,  the  solution  proceeds  with  two  unknown  functions  including 
the  potential,  q{z),  and  f{z,k)  which  is  a  solution  to  a  one-dimensional  Schrodinger 
equation  of  the  form 

/"  +  {k^  -  <l)f  =  0.  (1.1) 

The  above  equation  is  written  and  solved  in  terms  of  the  function  /  alone  by  using 
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the  following  “trace  formula”  to  eliminate  q{z)  in  equation  (1.1), 

2z 

q{z)  =  -  kR{k)f\z,k)dk  (1.2) 

TT  J-ca 

where  the  reflection  coefficient,  R{k),  is  the  input  data  and  k  is  the  wavenumber. 

The  trace  method  was  first  applied  to  the  shallow-water  ocean  acoustic  problem 
by  Stickler  and  Deift[77].  They  showed  how  to  recover  the  sound  speed  in  a  stratified 
ocean  where  the  frequency  was  chosen  so  that  no  proper  modes  were  excited.  As 
input  data  to  the  inversion  they  required  measurements  of  the  normal  derivative 
of  a  harmonic  point  source  pressure  field  at  the  pressure-release  surface.  Stickler[78] 
removed  the  requirement  of  measuring  the  normal  derivative  of  pressure,  while  keeping 
the  restriction  of  no  trapped  modes,  to  recover  the  sound  speed  and  density  profiles 
in  the  sediment  using  trace  methods.  He  used  measurements  of  the  sound  field  versus 
range  at  fixed  depth  for  two  distinct  frequencies  as  input  data. 

In  contrast  to  the  trace  methods  which  use  the  reflection  coefficient  directly,  the 
Gelfand-Levitan  (G-L)  approach  uses  the  Fourier  transform  of  the  reflection  coefficient 
as  input  data  to  a  Fredholm  integral  equation  of  the  second  kind.  Merab[50]  studied 
several  different  solution  methods  to  the  G-L  integral  equation  for  a  deep  water 
experiment.  The  reflection  coefficient  is  obtained  by  measuring  both  the  direct  and 
bottom  reflected  field  for  a  single-frequency  point  source  in  deep  water  and  applying 
the  Hankel  Transform  inversion  technique  described  in  [25].  The  Fourier  transform 
of  the  measured  reflection  coefficient  is  then  used  as  the  input  for  the  following  G-L 
type  integral  equation: 

K{z,y)  +  R{z +  y)  +  [  R{z' +  y)K{z,  z')dz' =  0,  y  <  z,  (1.3) 

where  R{z)  is  the  Fourier  transform  of  the  reflection  coefficient  R{k),  and  the  kernel, 
K{z,z)^  is  sought.  Given  K{z,z),  there  exists  a  simple  relation  to  find  the  desired 
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potential  q{z). 


(1.4) 


q{z)  =  2 


dK{z,z) 

dz 


Merab  illustrates  the  application  of  this  approach  using  synthetic  data  for  a  known 
geoacoustic  environment,  taking  care  to  avoid  having  trapped  modes  in  the  sediment. 

Each  of  the  above  methods  allows  for  the  presence  of  bound  states  or  trapped 
modes  in  the  solution  of  the  inverse  problem.  However,  in  each  of  the  cited  examples, 
cases  were  studied  where  bound  states  were  not  present.  This  was  done  because  the 
bound  states  further  complicate  the  mathematics  and  also  require  additional  mea¬ 
surements  to  be  made.  Further,  in  the  shallow- water  environment  of  interest,  it  is 
often  very  difficult  to  measure  the  plane- wave  reflection  coefficient  due  to  the  presence 
of  multiple  reflections  between  the  surface  and  the  bottom.  It  is  believed  that  this 
problem  is  compounded  by  the  difficulty  in  making  accurate  measures  of  the  phase 
of  the  complex  pressure  experimentally  [64].  In  the  Modal  Mapping  Experiments 
(MOMAX)  described  in  chapter  4,  it  will  be  demonstrated  that  much  progress  has 
been  made  in  this  regard  and  the  determination  of  the  reflection  coefficient  for  real 
data  will  be  revisited.  Recently,  McLaughlin  and  Wang  [51]  have  suggested  a  way 
of  applying  trace  methods  to  obtain  bottom  sediment  properties  from  shallow  water 
data.  In  this  formulation,  trapped  modes  propagating  in  the  water  column  are  treated 
as  bound  states  in  the  inverse  problem  and  the  potential  includes  the  sound  speed  in 
the  water  column.  McLaughlin  and  Wang  use  the  depth-dependent  Green’s  function 
obtained  by  Hankel  transforming  measurements  of  the  complex-pressure  field  as  their 
starting  point.  A  straightforward  relation  between  the  Green’s  function  and  the  re¬ 
flection  coefficient  is  applied  to  obtain  the  input  data  for  the  inversion.  To  account 
for  the  discrete  spectrum  present  in  the  Green’s  function  as  it  applies  to  the  trace 
formula,  McLaughlin  applies  a  Darboux  transform  [17]  to  systematically  remove  the 
bound  states  from  the  problem.  The  result  of  each  transform  step  yields  a  modified 
eigenvalue  equation  with  one  less  bound  state.  All  other  bound  states  remain  un¬ 
changed,  and  a  new  potential  results.  The  Darboux  transform  is  applied  iteratively 
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for  each  of  the  bound  states  and  information  is  stored  at  each  step.  After  removing 
all  the  bound  states  from  the  original  eigenvalue  problem,  the  inversion  proceeds  as 
if  no  bound  states  were  present.  A  modified  reflection  coefficient  is  found  for  the  now 
continuous  spectrum  and  used  as  input  data  to  the  trace  formulas  to  solve  for  an 
intermediate  potential.  The  full  potential  is  obtained  from  the  intermediate  result  by 
reversing  the  Darboux  transformation  and  adding  the  bound  state  information  back 
in.  Although  not  considered  further  in  this  thesis,  the  Darboux  transform  appears 
of  interest  for  removing  the  effects  of  trapped  modes  in  the  sediments  avoided  by 
Merab.  Further  discussion  of  the  Darboux  transform  for  potential  application  to  this 
problem  can  be  found  if  Deift  and  Trubowitz  [17]  as  well  as  Beals  et  al.  [6]. 

1.4  Experimental  Approach 

In  the  next  section,  the  concept  of  wavenumber  evolution  is  illustrated  using  a  nu¬ 
merical  example.  Results  show  the  dependence  of  horizontal  wavenumber  on  local 
waveguide  properties.  Following  the  numerical  results,  the  experimental  method  for 
extracting  local  wavenumber  content  is  briefly  introduced. 

1.4.1  Numerical  Results 

As  stated,  the  input  data  for  both  the  perturbative  inverse  method  and  Gelfand- 
Levitan  method  is  the  spectral  data  for  the  point  source  pressure  field.  The  under¬ 
lying  premise  is  to  measure  the  spatial  variability  of  the  wavenumber  spectrum  for 
a  point  source  in  a  stratified  medium  to  recover  the  geoacoustic  properties  of  the 
seabed  in  three  dimensions.  As  part  of  this  work,  both  at-sea  and  numerical  exper¬ 
iments  have  been  conducted  to  test  these  concepts.  Much  of  the  work  was  driven 
by  the  notion  that  perturbative  inversion  methods  would  be  used  to  infer  sediment 
properties  using  eigenvalue  information  only.  Thus,  methods  were  investigated  for 
extracting  and  mapping  the  wavenumber  spectrum  of  the  normal  mode  field  for  com- 
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Figure  1-2:  Three-dimensional  environmental  model.  Sound  speed  profile  is  shown  in 
upper  left.  Cross-range  and  down-range  bathymetry  profiles  are  shown  in  right  upper 
plots.  Full  3-D  model  is  shown  in  bottom  plot  with  source  located  at  30  meters  depth 
in  center  of  model.  Sediment  sound  speed  profile  tracks  the  bathymetry. 

plex  shallow-water  waveguide  environments  whose  acoustic  properties  vary  in  three 
spatial  dimensions[7].  For  the  waveguide  depicted  in  figure  1-2,  50  Hz  synthetic  pres¬ 
sure  data  was  generated  using  the  KRAKEN  normal  mode  program  [59]  as  shown 
in  figure  1-3.  Short-aperture  sliding  window  transforms  were  applied  along  radials 
to  obtain  estimates  of  the  Green’s  function  g{kr,r,6).  Estimates  where  made  using 
a  classical  spectral  estimator  based  on  the  FFT  with  a  Hanning  window  applied  to 
the  data  for  each  aperture.  Wavenumbers  corresponding  to  the  first  two  modes  were 
extracted  and  plotted  as  a  function  of  range  and  theta  to  create  the  modal  maps  in 
Figure  1-4.  These  maps  represent  the  ideal  case  w'here  modal  information  is  obtained 
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Pressure  Field  for  Three-Dimensional  Waveguide 


Figure  1-3:  Plan  view  of  50  Hz  Nx2D  acoustic  pressure  field  generated  nsing  KRAKEN 
for  3-D  environmental  model. 
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on  a  fully  populated  spatial  array.  The  wavenumber  estimates  obtained  represent 
the  true  wavenumber  values  well,  but  with  low  spatial  resolution  due  to  the  length 
of  the  aperture  used  at  each  range  step.  To  improve  the  spatial  resolution  of  the 
spectral  measurements,  a  portion  of  this  thesis  investigates  a  high-resolution  spectral 
estimation  technique  to  obtain  modal  wavenumber  information  using  short  aperture 
data  [9].  This  will  be  discussed  in  detail  in  chapter  3. 

1.4.2  Modal  Mapping  Experiments 

To  test  these  methods  on  real  data,  the  Modal  Mapping  Experiments  (MOMAX) 
were  conducted  in  March  1997  and  again  in  October  2000  off  the  New  Jersey  coast 
[29]  along  with  MOMAX  II  held  in  February  1999  in  the  Gulf  of  Mexico  [47].  In  these 
experiments,  drifting  buoys  equipped  with  a  single  hydrophone,  GPS  navigation,  and 
radio  telemetry  were  used  to  measure  the  acoustic  field  from  a  point  source  emitting 
several  discrete  tones  in  the  frequency  range  20-475  Hz.  High-resolution  measure¬ 
ments  of  the  sound  field  using  several  buoys  were  used  to  create  two-dimensional 
synthetic  aperture  horizontal  arrays.  Modal  maps  were  generated  for  both  single  re¬ 
ceiver  arrays,  and  for  arrays  formed  by  merging  data  from  two  calibrated  receiving 
systems[8]  [10].  Figure  1-5  shows  one  result  from  MOMAX  comparing  the  measured 
field  to  a  synthetic  field  generated  using  the  recovered  sound  speed  profile  obtained 
using  a  perturbative  inversion  algorithm  [62]  as  input  data  for  an  acoustic  propa¬ 
gation  model.  These  results  illustrate  the  utility  of  the  measured  Green’s  function 
spectrum  for  characterizing  the  spatial  variability  of  the  acoustic  waveguide  and  in¬ 
ferring  its  geoacoustic  properties.  Further,  the  phase  of  the  measured  acoustic  data 
from  MOMAX  appears  to  be  very  robust  and  shows  promise  for  extracting  the  re¬ 
flection  coefficient  from  the  depth-dependent  Green’s  function.  For  the  case  of  no 
trapped  modes  in  the  sediment,  these  measurements  may  be  used  for  inversion  using 
the  Gelfand-Levitan  method. 
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Figure  1-5:  Model  data  comparison  for  MOMAX  experiment  at  75  Hz 

1.5  Outline 

Chapter  2  presents  an  overview  of  the  theory  of  wave  propagation  in  a  horizontally 
stratified  shallow-water  waveguide.  The  point-source  solution  is  given  in  terms  of  the 
depth-dependent  Green’s  function.  The  relationship  between  the  Green’s  function 
and  plane-wave  reflection  coefficient  of  the  bottom  is  given  highlighting  the  Green’s 
functions  dependence  on  bottom  properties.  Doppler  effects  on  the  Green’s  func¬ 
tion  due  to  source/receiver  motions  is  also  examined.  The  final  part  of  the  chapter 
describes  the  inverse  problem.  The  perturbative  inverse  method  is  reviewed  as  appli¬ 
cable  to  the  MOMAX  experiments.  Finally,  the  Gelfand-  Levitan  inverse  method  is 
reviewed  with  emphasis  placed  on  applications  using  real  data. 

Chapter  3  presents  data  and  results  using  synthetic  data  generated  for  range- 
dependent  shallow-water  environments.  Data  was  provided  as  part  of  an  ONR/SPAWAR 
sponsored  inversion  techniques  workshop  [40]  that  occurred  in  May,  2001.  A  high- 
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resolution  wavenumber  estimator  to  extract  range-dependent  spectra  was  implemented 
for  this  work.  In  this  chapter,  the  estimator  will  be  described  and  its  statistical  char¬ 
acteristics  examined  for  application  to  the  MOMAX  data.  A  related  study  will  be 
described  that  examines  the  effect  of  sound  speed  perturbations  in  the  water  column 
on  wavenumber  measurements. 

Chapter  4  describes  in  detail  the  modal  mapping  experiments.  An  overview  of  the 
experimental  configuration  is  given  along  with  a  description  of  the  data  processing 
scheme.  An  overview  of  the  operational  environment  for  each  experiment  is  then 
given  along  with  supporting  measurements  of  the  water  column  and  bathymetry. 
Finally,  some  example  measurements  and  preliminary  analysis  of  the  MOMAX  data 
are  discussed. 

Chapter  5  presents  analysis  of  the  MOMAX  experimental  data.  Modal  evolution 
is  examined  for  the  tracks  created  by  the  drifting  buoys.  Analysis  is  presented  for 
the  case  of  a  single  propagating  mode  where  modal  evolution  can  be  determined 
analytically.  For  the  case  of  a  moving  source,  a  method  for  measuring  the  group 
velocity  is  presented.  MOMAX  data  are  inverted  using  the  perturbative  method  to 
obtain  compressional  wave  speed  profiles  in  the  sediment.  Finally,  applicability  of  the 
MOMAX  data  to  the  Gelfand-Levitan  inversion  method  is  discussed  with  suggestions 
for  further  work. 

The  final  chapter  summarizes  the  present  work  with  emphasis  placed  on  new 
contributions.  Suggestions  for  future  efforts  conclude  the  chapter. 
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Chapter  2 


Review  of  Sound  Propagation  in  a 
Shallow- Water  Waveguide 


2.1  Forward  Problem 


The  governing  equations  for  sound  propagation  in  shallow  water  are  reviewed  for  the 
case  of  a  horizontally  stratified  medium.  The  shallow  water  model,  as  shown  in  figure 
1-1,  is  comprised  of  a  stratified  fluid  layer,  with  sound  velocity  and  density  both 
functions  of  depth,  overlying  a  seabed  located  at  depth  zi,  made  up  of  multiple  layers 
whose  properties  vary  with  depth. 

The  starting  point  is  the  full  time-dependent  acoustic  wave  equation  in  three- 
dimensions  for  a  time-harmonic  source  with  radial  frequency  w,  written  [28], 


Po(r)V 


1 

Po(r) 


VP(r,t) 


1  d^P(r,t) 
c^(r)  df^ 


—A'k5{v  —  To)e 


(2.1) 


where  r  represents  the  spatial  coordinates  {x,y,z)^  Fq  is  the  source  position,  t  is 
time,  5(r)  is  the  Dirac  delta  function,  and  V  =  -I-  The  density  po  and 

sound  speed  c  are  functions  of  position  and  suggest  the  influence  of  the  propagation 
environment  on  the  solution  to  the  wave  equation.  Using  the  following  definition  for 
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the  Fourier  transform  in  going  from  the  time  domain  to  the  frequency  domain, 

T.r.{x{t)}  =  X{u)  =  ^  /  x{t)e^^dt,  (2.2) 

v27r  J-oo 

and  applying  it  to  (2.1),  the  time-independent  Helmholtz  equation  results, 

Pc.(r)V-  -^^Vp(r)  +A;2(r)p(r)  = -47r(5(r-ro),  (2.3) 

Here,  angular  frequency  is  defined  by  w  =  27r/,  for  frequency  /,  and  A:(r)  =  oj/c{v)  is 
the  total  acoustic  wavenumber. 

The  wave  equation  can  be  simplified  for  our  case  by  considering  the  solutions 
in  regions  of  constant  density.  For  a  stratified  waveguide  where  density  is  constant 
within  each  layer,  equation  (2.3)  takes  the  form 

[V2  -h  A:2(r)]  p(r)  =  -AT,5{x)5{y)5{z  -  z,),  (2.4) 

for  the  source  located  at  the  origin  of  the  (x,y)-plane  and  at  a  depth  Zg. 

2.1.1  Point-source  spectrum  and  plane-wave  reflection  coef¬ 
ficient 

For  a  waveguide  with  properties  that  vary  only  as  a  function  of  depth,  the  solu¬ 
tion  to  (2.4)  proceeds  by  reducing  the  equation  to  one  in  the  depth  coordinate  only. 
Appropriate  boundary  conditions  are  then  applied  at  the  air-sea  and  water-bottom 
interfaces.  Defining  the  two-dimensional  inverse  Fourier  transform  operator  in  spatial 
coordinates  by  [28], 

I.T.r.{f{x,y)}  =  F{k^,ky)  =  —  /  /(x,  (2.5) 

ZTT  J —oo  J —oo 
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and  applying  it  to  both  sides  of  (2.4),  the  following  equation  results, 


2S(^Z  ■2'o). 


(2.6) 


Here  is  the  vertical  wavenumber  defined  as  kz  =  ~  kx  and  ky  the 

horizontal  wavenumbers  with  respect  to  a  Cartesian  coordinate  system  in  wavenumber 
space.  In  the  above,  it  is  required  that  lTn{kz}  >  0  in  order  to  ensure  the  field  is  finite 
as  z  — >  oo.  The  solution  g{kx,  ky]  z,  Zo)  is  the  depth-dependent  Green’s  function  or 
point-source  spectrum  for  the  depth-separated  wave  equation.  The  Green’s  function 
is  related  to  the  acoustic  pressure  field  through  the  conjugate  Fourier  transform  pairs 
given  by: 


p(x,  p-,z,zj  =  ^  r  r  g(k„  k,-,z,  (2.7) 

ZTT  J —oo  J  —  CO 

g{kx,  ky]  z,Zo)=  -^  [  j  p(x,  y]  z,  (2.8) 

ZTT  J  —oo  J —oo 

The  solution  to  the  depth-dependent  equation  (2.6)  is  obtained  by  incorporating 
into  the  total  solution,  the  independent  solutions  urikx,  ky]  z)  and  UB{kx,  ky]  z)  which 
satisfy  the  top  and  bottom  boundary  conditions,  respectively.  The  general  form  of 
the  depth-dependent  Green’s  function  satisfying  both  boundary  conditions  can  be 
written. 


g(k„k,-,z.z,)  =  _'^r(kz.ky-,z^UB{)=„k,-,z>) ,  „  <  ^ 


(2.9) 


where  2<  =  min{z,Zo),  =  max{z,Zo)  and  W(zo)  is  the  Wronskian,  given  in  terms 
of  the  solutions  ut  and  ub  evaluated  at  the  source  location  by. 


W(zo)  =  Ut{Zo)u'q{Zo)  -  u'.j^{zo)ub{zo)-  (2.10) 


For  an  isovelocity  water  column,  the  solutions  ut  and  ub  can  be  expressed  in 
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terms  of  up-going  and  down-going  plane  waves  along  with  their  boundary  interactions 
characterized  by  the  plane-wave  reflection  coefficients  RriK)  and  Rb{K)  at  the  top 
and  bottom  interfaces.  Expressing  the  reflection  coefficients  as  a  function  of  horizontal 
wavenumber,  kr,  the  two  independent  solutions  are  given  by, 

UT{kr\ z)  =  A  +  RT{K)e’^^^] ,  (2.11) 

UB{kr,  z)  =  B  -  RB{kr)e^^^] ,  (2.12) 

for  which  the  Wronskian  becomes, 

W  =  2ik,AB(l  -  (2.13) 


where  A  and  B  are  arbitrary  constants  and  kr  =  Substituting  these  expres¬ 

sions  into  (2.9),  the  depth-dependent  Green’s  function  for  a  shallow  water  waveguide 
is 


g{kr\ z,  Zo)  = 


+  RB{kr)e^^'^^^»  -h  RT{kr)e 

—  RT{kr)RB(k.r)e^^^^^>>\ 


~ikz\z~Zo\ 


(2.14) 

This  expression  for  the  Green’s  function  is  a  spectral  representation  of  the  propagating 
field  in  terms  of  horizontal  wavenumber  and  shows  explicitly  its  dependence  on  the 
boundary  conditions  of  the  waveguide.  For  the  ocean  acoustic  problem,  the  air-sea 
interface  is  approximated  by  a  pressure-release  surface  with  boundary  condition  given 
by  Rt  =  —  1.  The  depth-dependent  Green’s  function  is  then  an  explicit  function  of 
the  plane-wave  reflection  coefficient  at  the  water-bottom  interface.  For  example,  for 
a  stratified  seabed  with  multiple  layers,  as  illustrated  in  figure  (2-1),  Rb{k.r)  can  be 
written  as  a  sum  including  the  plane-wave  reflection  coefficients  calculated  between 
adjacent  layers  [44],  and  a  phase  term  that  depends  on  the  layer  depths. 
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Figure  2-1:  Layered  sediment  model  used  for  calculation  of  plane-wave  reflection 
coefficient. 


Rb{kr)  Boi 

+ 

+ 

1 

+ 

+ 

Ri2 

1 

+ 

+ 

D-1  ^-2iA:n2{dAf 

Rn{n+\)^ 

(2.15) 


where  Rij  is  the  reflection  coefficient  between  two  adjacent  layers  given  by 

(2.16) 

Pjr^zi  -f"  Pi^zj 

From  (2.15)  it  is  clear  that  Rb  depends  solely  on  the  sediment  properties  of  the 
bottom.  Consequently,  measurements  of  the  bottom  reflection  coefficient  are  often 
used  as  a  basis  for  the  geoacoustic  inverse  problem  such  as  in  the  trace  methods  and 
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the  method  of  Gelfand-Levitan  previously  discussed.  In  the  stratified  shallow-water 
case,  the  relationship  between  Ri,{kr)  and  g{kr)  given  by  (2.14)  shows  the  explicit 
dependence  of  the  spectrum  on  the  bottom  properties.  This  relationship  is  the  basis 
for  geoacoustic  inverse  problems  using  point  source  acoustic  field  data  in  shallow- 
water.  Of  particular  interest  for  the  inverse  problem  is  the  discrete  portion  of  the 
wavenumber  spectrum  g  given  by  (2.14)  when  the  denominator  goes  to  zero.  Making 
the  substitution  Rt  =  —1, 

1  +  =  0.  (2.17) 

It  will  be  shown  that  this  equation  takes  the  exact  form  of  the  characteristic  equation 
for  the  inhomogeneous  depth-dependent  wave  equation  for  a  waveguide  with  pressure- 
release  top  and  horizontally  stratified  bottom  [28].  Solutions  to  (2.17)  can  be  used  as 
input  data  to  the  perturbative  inverse  methods. 


2.1.2  Cylindrical  coordinates:  The  Hankel  Transform 

A  more  natural  representation  of  (2.7)  for  a  point  source  in  a  stratified  waveguide  is 
obtained  by  transforming  the  problem  to  cylindrical  coordinates.  For  an  axisymmetric 
problem,  transformation  to  cylindrical  coordinates  has  the  effect  of  reducing  the  two- 
dimensional  conjugate  Fourier  transform  pair  relationship  of  (2.7)  and  (2.8)  to  the 
one-dimensional  Hankel  transform  pair  relationship  with  conjugate  variables  kr  and  r. 
The  Cartesian/cylindrical  transform  relations  in  the  spatial  and  wavenumber  domains 
are  given  by  the  following  identities. 


k^  —  kf  cos  Q^, 
ky  =  kr  sin  a, 

-  yjkl  +  kl. 


X  =  r  cos  P, 
y  =  r  sin  /?, 
r  =  -f  y'^, 


(2.18) 


where  the  angles  oc  and  P  are  as  shown  in  figure  (2-2). 


Using  these  transforms. 
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Figure  2-2:  Angle  definitions  for  coordinate  transforms  in  Cartesian  and  wavenumber 
space 

equation  (2.7)  becomes 

1  roo  r2TT 

p(nz,z„)  =  Tr-/  /  (2.19) 

2jT^  «/ 0  « 0 

Using  the  integral  representation  for  the  zero-order  Bessel  function  Jo, 

J^{z)  =  ^  (2.20) 

ZTT  JO 

equation  (2.19)  becomes  a  zero-order  Hankel  transform  which  makes  up  one  half  of 
the  conjugate  transform  pair  given  by, 


p(r;  Zo)  = 

roo 

1  ^ (^r ) 5  j 

JO 

(2.21) 

g{kr\ Z,  Zo) 

roo 

=  /  p{r;z,Zo)Jo{krr)rdr, 

Jo 

(2.22) 

Using  the  following  identity  relating  the  Bessel  function  to  the  Hankel  function, 

=  i  [//<'>((c,r)  +  Hf>(A:,r)]  . 
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the  one-sided  integral  of  (2.21)  can  be  written  as  the  two-sided  integral  [28] 


1 

p{r\ z,  Zo)  =  -  g{kr;  z,  Zo)Hy\krr)krdkr. 

Z  7—00 


(2.23) 


This  form  is  particularly  convenient  when  approximating  the  Hankel  function  by  its 
asymptotic  form  for  k^r  1  for  which  the  above  integral  and  its  inverse  take  the 
form  of  Fourier  transforms.  Using 


in  (2.23)  gives, 


with  the  inverse, 


g  i7r/4  ^oo  j — 

P{T 'i  '^5  ^  y™ .  /  k^E  dk^pj 

v27rr  7-00 


(2.24) 


g^^/4  rOO 

g{kr;z,Zo)  ~  J^^p{r;z,Zo)Vre~^’'^^dr. 


(2.25) 


These  expressions  are  easily  incorporated  for  use  on  measured  acoustic  field  data 
where  the  wavenumber  spectrum  can  be  extracted  using  methods  based  on  the  Fast 
Fourier  Transform  (FFT)  or  other  computational  techniques  commonly  used  in  signal 
processing. 

Equations  (2.24)  and  (2.25)  are  often  used  as  the  basis  for  waveguide  characteri¬ 
zation  and  geoacoustic  inversion  techniques  based  on  spatial  measurements  of  point 
source  acoustic  wave  fields. 
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2.1.3  Normal  Mode  Representation 

An  alternative  solution  to  the  standard  wave  equation  given  by  (2.4),  and  written  in 
cylindrical  coordinates  as  [28] 

[V^  +  P(r)]  p(r)  =  -AJ-^5{e)6{z  -  z^),  (2.26) 

is  obtained  through  the  assumption  of  a  separable  solution  of  the  form, 

=  '^an{zo)'^n{z)Rn{r),  (2.27) 

n 

where  r  =  {r,0,z),  i?„(r)  and  ’5'n('2)  are  radial  and  depth  eigenfunctions  and  an{zo) 
is  amplitude.  The  important  result  here  is  that  for  constant-density  layers,  the 
satisfies  the  eigenvalue  equation 

+  A:L»n(z)  =  0;  =  kHz)  -  kl  (2,28) 

along  with  its  boundary  conditions,  and  kn  is  the  eigenvalue  for  the  mode  n.  For  a 
homogeneous  fluid  layer  bounded  from  above  and  from  below  by  horizontally  strat¬ 
ified  media  with  plane-wave  reflection  coefficients  Rt  and  Rb,  respectively,  the  ho¬ 
mogeneous  solution  to  the  depth  equation  can  be  written  in  terms  of  up-going  and 
down-going  plane  waves. 

^(z)  =  ^rf(z)  -h  ^„(^)  =  (2.29) 


where  A  and  B  are  arbitrary  constants  determined  by  the  boundary  conditions.  Using 
this  expression,  the  reflection  coefficients  at  the  top  and  bottom  surfaces  can  be 
written  as 


Rx  — 


'^diz)  I  _  A 


(2.30) 
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(2.31) 


Rb  — 


I  _  _  ^p-2ik:,Zi 

'l'd(.z) A 


These  two  equations  can  then  be  combined  to  yield  the  characteristic  equation  for 
the  depth  dependent  equation  given  by 


1  -  =  0. 


(2.32) 


Recalling  equation  (2.14),  this  is  exactly  the  form  of  the  denominator  for  the  shallow- 
water  depth-dependent  Green’s  function  given  in  equation  (2.17)  for  Rt  =  —1.  Poles 
of  the  Green’s  function  thus  correspond  exactly  to  eigenvalues  for  the  depth-dependent 
eigenvalue  equation. 

Following  the  treatment  of  Frisk  [28],  the  radial  component  of  the  solution  (2.27) 
can  be  shown  to  be  a  solution  to  Bessel’s  equation  with  a  line  source  driving  term 
and  has  the  form, 

Rn{r)  =i7rH^^'>{knr).  (2.33) 

The  modal  representation  of  the  acoustic  pressure  field  can  then  be  written, 

CX) 

p{r,  z)  =  in  an{zo)'^niz)H^^\knr).  (2.34) 

n=l 

Finally,  solving  for  an{z)  for  z  =  Zo  using  the  closure  relation  for  eigenfunctions  of  a 
proper  Sturm-Liouville  system  as. 


(^n(Zo)  ~  ^n('^o)/p('2o)) 


and  the  modal  representation  becomes, 


PK^o)  n=l 


(2.35) 


(2.36) 
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In  the  far  field,  the  asymptotic  form  of  the  Hankel  function  can  be  used  to  obtain. 


p(r,  z)  - 


p{Zo) 


00  piknT 


(2.37) 


2.1.4  Range-Dependent  Modal  Representation 

In  the  previous  development  of  the  modal  solution  to  the  depth-dependent  problem, 
a  range-independent  medium  was  assumed.  Pierce  [57]  first  introduced  the  range- 
dependent  formulation  for  normal  modes  where  he  assumed  small  variations  in  the 
propagation  medium  with  range  and  that  individual  modes  do  not  transfer  energy 
between  one  another.  Thus,  by  representing  the  solution  at  any  given  range  as  a  sum 
of  local  modes,  where  the  modes  satisfy  the  local  boundary  conditions,  the  adiabatic 
mode  formula  can  be  derived  and  is  given  by. 


p{r,  z) 


iv/i 

p{Zo) 


oo 


z) 


k„(r')dr' 
y/o  kn{r')dr' 


(2.38) 


where  the  integral  over  dr  in  the  denominator  of  the  mode  sum  has  been  included  to 
ensure  reciprocity  between  source  and  receiver.  Through  this  formulation  and  appli¬ 
cation  of  the  transform  relations  between  the  depth-dependent  Green’s  function  and 
the  complex-pressure  field,  range-dependent  spectra  can  be  extracted  by  performing 
the  wavenumber  integral  over  local  apertures  of  finite  length. 

It  should  be  noted  that  although  adiabatic  mode  theory  is  a  convenient  way  to  rep¬ 
resent  the  range-varying  nature  of  the  local  modes  for  a  spatially  varying  waveguide, 
it  is  not  a  necessary  condition  for  the  extraction  of  range-dependent  spectra.  It  will 
be  demonstrated  that  effective  wavenumber  estimates  can  be  obtained  for  regions  of 
locally  range-independent  environments  bounded  by  regions  of  the  environment  where 
adiabatic  mode  theory  does  not  apply  and  mode  coupling  occurs.  Similarly,  for  the 
case  of  weak  internal  waves  in  a  range-independent  waveguide,  effective  estimates 
of  the  propagating  modes  can  be  made  where  mode  coupling  occurs.  In  this  case 
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it  will  be  shown  that  the  higher  order  modes  are  relatively  stable  indicating  range- 
independent  bottom  properties  and  that  energy  transfered  amongst  the  lower-order 
modes  is  due  to  range-dependence  in  the  water  column. 


2.1.5  Effects  Due  to  Source/Receiver  Motion 

In  the  MOMAX  experiments  to  be  described,  the  source/receiver  geometry  was  de¬ 
pendent  on  the  drift  tracks  of  the  individual  buoys  and/or  whether  the  source  was 
moored  or  being  towed.  As  such,  any  effect  these  configurations  might  have  on  the 
measured  wavenumber  spectra  should  be  addressed.  In  the  typical  range  separation 
experiment,  it  has  been  recognized  that  towed  source  motion  imparts  a  shift  in  the 
measured  wavenumber  spectrum.  Rajan  et.  al.  [67]  corrected  for  this  effect  in  mea¬ 
surements  made  in  the  Hudson  Canyon  area  of  the  Atlantic  Ocean,  where  a  source 
was  towed  along  a  radial  at  1.5  m/s.  Hawker  [38]  and  Schmidt  and  Kuperman  [73] 
give  detailed  analyses  of  the  moving  source  problem  and  derive  independent  expres¬ 
sions  for  the  expected  Doppler  shifts.  Each  of  their  methods  will  be  reviewed  here 
and  applied  to  the  MOMAX  measurements  in  a  following  chapter. 

Hawker  treats  the  case  for  a  source  moving  along  a  constant  trajectory  with  con¬ 
stant  speed  past  a  stationary  receiving  array.  For  the  geometry  in  Figure  2-3  he 
derives  a  modal  solution  to  the  wave  equation  that  includes  the  geometry  and  source 
motion.  Through  his  analysis,  he  arrives  at  a  modal  solution  which  depends  on  source 
velocity  and  angle  to  receiver  along  with  modal  group  velocity. 


p{r,  z) 


- — - —  \  - — exp 


KR 


V 


-sine) 

Vn  J  i 


(2.39) 


where  R  and  6  are  the  slant  range  and  angle  between  source  and  receiver,  k°  is 
the  modal  eigenvalue  for  no  source  motion,  and  is  the  group  velocity  of  mode  n. 
Comparing  the  result  of  Hawker  given  by  (2.39),  to  the  mode  sum  for  a  stationary 
source  given  by  (2.37),  the  horizontal  wavenumbers  in  the  moving  source  problem  are 
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Figure  2-3:  Source/Receiver  geometry  considered  by  Hawker  for  moving  source  prob¬ 
lem  [38]. 

shifted  by  an  amount  equal  to  the  ratio  of  the  source  speed  to  modal  group  velocity. 
For  the  geometry  in  figure  2-3,  the  shifted  wavenumbers  are  given  by, 

K  -  •  (2.40) 

A  later  paper  by  Schmidt  and  Kuperman  [73]  treats  a  similar  problem,  but  allows 
for  more  general  source/receiver  geometry.  In  their  analysis,  the  source  and  receiver 
are  constrained  to  the  same  horizontal  plane,  but  can  move  with  constant  velocity  at 
arbitrary  angles  with  respect  to  one  another  as  illustrated  in  figure  2-4.  Under  these 
conditions,  and  assuming  a  spatially  independent  horizontally  stratified  medium,  the 
spectral  and  modal  representations  of  the  pressure  field  are  derived.  The  Green’s 
function  is  shown  to  be  a  solution  to  a  depth-dependent  wave  equation  as  in  the 
stationary  case  of  (2.6)  but  evaluated  at  a  Doppler  shifted  frequency  given  by 

ui*  =  UJ  +  K-  Vs,  (2.41) 
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Figure  2-4:  Geometry  used  by  Schmidt  and  Kuperman  for  analysis  of  Doppler  shifts 
due  to  source/receiver  motions  [73]. 

where  Vg  is  the  source  velocity.  The  time-domain  solution  to  the  wave  equation  for 
a  moving  source  is  then  shown  to  be  given  in  terms  of  a  depth-dependent  Green’s 
function  that  depends  on  the  source  motion.  The  receiver  motion  is  introduced  inde¬ 
pendently  and  is  shown  to  introduce  a  frequency  shift  in  the  exponential  term  of  the 
time-domain  solution,  but  not  in  the  integration  kernel.  For  both  a  moving  source 
and  moving  receiver,  the  time-domain  solution  to  the  wave  equation  can  be  written 
as  [73] 

p(ro  +  Vrt,  z,  t)  ~  J  g(kr,  z-  uJo  +  K-  (2.42) 

where  the  receiver  is  moving  with  a  velocity  vector  Vp  given  by  r  =  Fq  4-  This 
result  of  Schmidt  and  Kuperman  shows  that  the  effects  of  moving  sources  and  re¬ 
ceivers  is  not  reciprocal.  Only  source  motion  causes  a  Doppler  frequency  shift  in  the 
propagating  modes  of  the  field  for  a  harmonic  source,  while  both  source  and  receiver 
motion  contribute  to  the  Doppler  wavenumber  shift  in  the  phase  term.  In  analysis 
of  the  experimental  data  for  this  work,  observations  of  the  modal  shift  predicted  by 
these  models  will  be  explored  and  a  method  for  extracting  the  group  velocity  directly 


from  measurements  considered. 


2.2  Inverse  Problem 

The  preceding  sections  describe  the  forward  problem  for  predicting  acoustic  propaga¬ 
tion  as  a  solution  to  the  depth-dependent  wave  equation  given  the  boundary  condi¬ 
tions.  The  inverse  problem  of  interest  is  to  determine  the  properties  of  the  bounding 
media  through  measurements  of  the  propagating  acoustic  field.  The  data  used  in 
the  inverse  methods  to  be  described  is  the  spectral  representation  of  the  propagat¬ 
ing  acoustic  field.  It  was  shown  in  equation  (2.14)  how  this  representation  contains 
information  about  the  bottom  properties  of  interest. 

2.2.1  Perturbative  Inversion 

The  perturbative  inverse  method  is  an  iterative  method  that  is  based  on  making 
small  perturbations  to  a  starting  background  model.  Here,  the  formalism  of  Rajan  et 
al.  [62]  is  followed  for  formulating  the  problem  as  a  Fredholm  integral  equation  that 
can  be  solved  using  linear  inverse  theory. 

The  starting  point  is  the  solution  to  the  depth-dependent  wave  equation  given 
a  predetermined  environmental  model.  The  normal  mode  solution  for  the  initial 
problem,  assuming  a  piecewise  constant  density  profile  is  written, 

{£;  +  (2.43) 

where, 

k^°'>{z)=u}/coiz)  (2.44) 

is  the  wavenumber  for  the  background  sound  speed  profile  Co{z),  'F[)]^(z)  is  the  mth 
normal  mode  function,  and  is  the  mth  mode  eigenvalue  for  the  background  prob¬ 
lem.  Given  this  model,  a  perturbation  in  the  background  sound  speed  model  by  Ac(z) 
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gives  the  perturbed  wavenumber, 


k{z)  =  oj/[co{z)  +  Ac(z)].  (2.45) 

This  perturbation  in  wavenumber  yields  a  corresponding  perturbation  to  the  eigen¬ 
values  that  when  linearized  can  be  expressed  by, 

Aft™  =  (2ftm)-'{»S>|ft?l»S>),  (2.46) 

where  brackets  denote  integration  over  depth,  and  H  is  the  wavenumber  perturbation 
squared  and  linearized. 

H  =  Ak^{z)  =  e{z) 

=  +  (2«) 
=  —  A:^°^^(2Ac/co)  for  ^cjcg  <<  1. 

Substituting  equation  (2.47)  into  (2.46)  and  expressing  the  depth  integral  indicated 
by  the  brackets  explicitly,  the  mode  number  perturbations  can  be  written  as  the 
integral  equation, 

r (2.48) 

km  ^o\^) 

In  this  form,  a  given  small  perturbation  in  the  background  sound  speed  profile  results 
in  a  small  change  in  the  modal  eigenvalues.  The  iterative  perturbative  inversion 
method  thus  makes  small  perturbations  to  the  background  sound  speed  profile  and 
compares  the  resulting  eigenvalues  to  the  measured  input  values.  The  background 
profile  is  updated  iteratively  until  the  eigenvalues  from  the  updated  profile  agree  with 
the  measured  input  data. 

Rajan  et.  al.  [62]  details  three  methods  for  solution  of  the  integral  equation  given 
by  (2.48)  from  which  the  regularization  method  is  shown  here.  This  method  applies  a 
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smoothness  constraint  on  the  solution  of  the  integral  equation  by  requiring  the  second 
derivative  of  the  sound  speed  to  be  continuous  across  layer  interfaces  in  the  sediment. 
The  smoothness  measure  thus  stated  is  defined  as, 


Using  the  above  constraint,  the  integral  equation  (2.48)  is  discretized  and  recast  as 
an  augmented  least-square  problem  with  the  objective  of  minimizing  a  cost  function 
given  by  the  squared  difference  of  as  given  below.  Input  data  for  the  problem 
is  the  difference  between  measured  eigenvalues  as  determined  by  taking  the  Hankel 
transform  of  complex  pressure  and  the  eigenvalues  calculated  for  the  background 
model.  The  data  is  a  vector,  d,  whose  length  is  determined  by  the  number  of  modes 
with  elements  given  by 

di  =  /\kjn  =  k,,  (2.50) 

where  km  are  the  measured  eigenvalues  and  ki  are  the  mode  eigenvalues  for  the  back¬ 
ground  problem.  The  unperturbed  values  in  the  integral  equation  (2.48)  can  be 
expressed  as  a  matrix,  G,  with  elements, 

=  (2-51) 

where  i  and  j  refer  to  mode  number  and  layer  depth  indices,  respectively.  The 
smoothness  constraint  is  then  expressed  by  the  double  difference  equation  given  by, 

S{Ac)  =  -  2Ac,  +  Ac,_i)"  =  Ac^HAc,  (2.52) 

j 
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with 


H  = 


1-21  00 
-2  5-4  1 

1-4  6-4 


0  0 


(2.53) 


0  0  0  0  1  -2  1 


In  the  above,  S'  is  a  vector  of  length  depending  on  the  number  of  depth  layers  the 
problem  has  been  discretized  with,  G  is  an  n  x  m  matrix  corresponding  to  the  number 
of  modes  and  number  of  layers.  The  solution  is  then  found  by  determining  the  Ac 
that  minimizes  the  following  equation  for  I  given  by 


I  =  AAc^HAc  +  (d  -  GAc)^(d  -  GAc),  (2.54) 


where  A  is  a  Lagrange  multiplier  representing  the  amount  of  smoothness  applied. 
Differentiating  (2.54)  with  respect  to  each  element  of  Ac,  the  problem  takes  the  form 

Ac  =  (G'^G  +  AH)-^G'^d.  (2.55) 

In  the  above,  when  A  =  0,  no  smoothness  constraints  are  applied  and  the  equation 
reduces  to  the  standard  least  square  solution  for  Ac  given  by 


Ac  =  G-^d.  (2.56) 

Thus  stated,  for  a  given  Akm,  the  minimization  problem  is  solved  to  determine  a 
perturbation  to  the  background  model  by  Ac.  The  cost  function  is  evaluated  for  the 
new  background  model  and  iteration  proceeds  until  agreement  is  obtained  between 
the  updated  model  and  the  input  data. 
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2.2.2  Gelfand-Levitan  Method 


The  Gelfand-Levitan  method  is  an  exact  method  for  geoacoustic  inversion  that  is 
based  on  the  analogous  representation  of  the  depth-dependent  acoustic  equation  with 
the  time-independent  Schrddinger  equation  of  quantum  scattering  theory.  The  time- 
independent  inhomogeneous  wave  equation  for  the  depth-dependent  variable  u{kr,z) 
can  be  written  as, 

+  k'^  -  k^^  u(kr,  z)  =  0.  (2.57) 

By  making  a  simple  wavenumber  translation  given  by  p?  —  k^  —  /c^,  and  q{z)  — 
kl  —  k'^iz),  the  equation  (2.57)  can  be  written, 

(^-^  + ~  q{z)^u{fi,z)  =  0,  (2.58) 


with  boundary  condition  0)  =  0.  This  is  easily  seen  to  take  the  form. 


which  is  the  time-independent  Schrodinger  equation  where  E  =  is  the  “energy” 
and  V{z)  =  q{z)  is  the  “potential”.  In  the  above  and  for  what  follows,  it  has  been 
assumed  that  the  density  is  held  constant  over  all  depths.  This  requirement  will  be 
relaxed  in  a  later  section. 

Considering  the  model  shown  in  figure  2-5,  and  following  closely  the  develop¬ 
ment  by  Chadan  and  Sabatier  [12],  the  scattering  solution  to  this  full-line  equation 
is  governed  by  the  asymptotic  boundary  conditions  for  2  — >  ±00: 

±i/iz 

u±{z)  ~  (2.60) 

V  27r  V  z'K 


dz'^ 


+  [E-V{z)]\u{E,z)  =  0, 


(2.59) 


U:^{z)  ~  T±(//) 


e±it^z 

7^^ 


z  — >  ±00. 


(2.61) 


The  two  independent  solutions  u±{z)  correspond  to  plane  waves  emanating  from 
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Figure  2-5:  1-Dimensional  inverse  scattering  problem  model.  Inverse  problem  is  to 
determine  the  potential  V(z)  of  the  scattering  region. 

2:  =  =FcxD  where  R:^  are  the  reflection  coefficients  in  the  directions  and  T±  are 
the  transmission  coefficients  in  the  ±z  directions.  The  reflection  and  transmission 
coefficients  make  up  the  full  scattering  matrix  of  the  Schrodinger  equation 

T.ifi)  ) 

Assuming  the  end  potentials  as  z  — >•  ±00  are  both  equal  to  zero,  the  full  scattering 
matrix  can  be  determined  from 

Thus,  the  geoacoustic  inverse  problem  becomes  one  of  determining  the  potential 
q{z)  from  the  reflection  coefficient  i?_  (//)  including  information  from  both  the  discrete 
and  continuous  spectrum.  Note  that  for  kr  >  k^o,  n  becomes  imaginary.  The  energy 
fj?  is  always  real  and  is  positive  for  //  real.  The  potential  q{z)  is  in  general  positive. 
Chadan  and  Sabatier [12]  continue  with  the  solution  of  the  one-dimensional  inverse 
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scattering  problem  making  the  following  assumption  about  the  potential: 


(1  +  \z\)V {z)dz  <  oo. 


(2.62) 


A  solution  to  the  inverse  problem  can  be  found  using  the  Fourier  transform  of  the 
reflection  coefficient, 

R{z)  -  —  /  (2.63) 

Zti  J —oo 

Where  R{ij,)  =  This  expression  is  inserted  into  the  Gelfand-Levitan  type 

equation, 


iF(z,  y)  +  J?(z  +  y)  +  f  R[z'  y)K[z,  z  )dz  —0,  y  z  (2.64) 

J  —OO 

which  can  be  solved  by  relatively  straightforward  numerical  methods.  It  should  be 
noted  that  Deift  and  Trubowitz  make  the  same  assumptions  about  the  potential  to 
solve  an  identical  inverse  problem  in  terms  of  their  trace  formula.  Merab  [50]  derives 
the  same  equation  for  unequal  end  potentials  at  z  =  ±co. 

The  theory,  as  stated,  is  quite  straightforward  and  requires  only  the  plane-wave  re¬ 
flection  coefficient  as  input  data  to  the  inversion.  However,  there  are  several  caveats 
that  must  be  addressed  before  application  of  the  Gelfand-Levitan  method  can  be 
applied  in  practice.  It  was  already  mentioned  that  the  development  assumed  the 
density  to  be  constant  in  the  problem.  This  assumption  allows  the  plane-wave  reflec¬ 
tion  coefficient  to  be  expressed  solely  as  a  function  of  sound  speed  ratios  at  the  layer 
interfaces.  However,  it  was  shown  in  equation  (2.16)  that  for  the  acoustic  problem, 
the  sound  speed  and  density  are  coupled  for  determining  the  reflection  coefficient. 
Merab  [50]  gives  a  derivation  for  the  case  of  continuous  variation  of  density  in  the 
sediment,  but  does  not  address  the  jump  discontinuity  at  the  water-bottom  interface. 
This  and  other  practical  issues  associated  with  solving  the  G-L  integral  equation  will 
be  discussed  along  with  some  numerical  examples  in  chapter  5  of  the  thesis.  A  goal  is 
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to  address  issues  relating  to  the  application  of  the  method  to  acoustic  pressure  data 
measured  at  sea. 
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Chapter  3 


Wavenumber  Estimation 
Techniques 

3.1  Introduction 

In  this  chapter  practical  aspects  of  extracting  the  discrete  horizontal  wavenumber 
spectrum  from  range-dependent  acoustic  data  are  considered.  Of  particular  impor¬ 
tance  are  problems  where  the  expected  local  eigenvalues  are  closely  spaced  in  the 
wavenumber  domain.  In  these  cases,  classical  spectral  analysis  methods  show  the 
limit  of  resolution  for  detection  of  two  closely  spaced  wavenumber  values  to  be  in¬ 
versely  proportional  to  the  data  length.  In  cases  where  data  are  limited  to  a  given 
aperture,  or  the  environment  changes  rapidly,  it  may  be  impossible  to  resolve  indi¬ 
vidual  modes.  This  is  demonstrated  in  figures  3-1-3-3  for  a  50  Hz  complex  pressure 
signal  measured  on  a  5  meter  grid  where  spectra  are  determined  for  apertures  of  1000 
m  and  2000  m  using  both  classical  and  high-resolution  methods.  More  details  about 
the  field  and  the  resulting  wavenumber  estimates  are  discussed  later  in  this  chapter. 
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tc3  (50  Hz) 


Range  (m) 

Figure  3-1:  50  Hz  synthetic  pressure  field  data  for  waveguide  with  unknown  range- 
dependent  sediment  properties.  1000  m  and  2000  m  data  segments  are  shown  as  used 
for  making  wavenumber  estimates. 
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Spectral  Estimates  (tc3  50  Hz  2000  m  Ap.) 


Figure  3-2;  Comparison  of  spectral  estimates  for  2000  m  data  aperture.  Classical 
PSD  result  is  solid,  high-resolution  result  is  dashed,  and  true  wavenumber  values  are 
dotted. 


Spectral  Estimates  (tc3  50  Hz  1000  m  Ap.) 


Figure  3-3:  Comparison  of  spectral  estimates  for  1000  m  data  aperture.  Classical 
PSD  result  is  solid,  high-resolution  result  is  dashed,  and  true  wavenumber  values  are 
dotted. 
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In  making  spectral  estimates  based  on  classical  methods,  Matlab’s  [60]  psd  func¬ 
tion  was  used  with  a  Hanning  window  applied  to  the  full  data  aperture  and  the 
number  of  FFT  points  specified  as  16384.  The  High-resolution  estimates  were  made 
using  the  modified-covariance  autoregressive  (AR)  estimator  of  Marple  [49]  with  AR 
model  orders  of  20  and  30  for  the  2000  m  and  1000  m  apertures,  respectively.  Having 
estimated  the  AR  model  coefficients,  the  spectra  were  determined  on  a  16384  point 
wavenumber  grid.  Full  details  of  the  AR  estimator  are  discussed  in  Marple  [49]  with 
a  summary  provided  below.  Before  proceeding,  it  should  be  emphasized  that  for  this 
work  the  data  to  be  extracted  from  the  wavenumber  spectral  plots  is  the  location  of 
the  wavenumber  peaks,  not  the  amplitudes.  Recognizing  that  the  variance  in  peak 
levels  for  AR  spectral  estimates  can  be  large  compared  to  other  methods  [4],  observed 
discrepancies  in  spectral  amplitudes  at  the  resolved  frequencies  between  the  classical 
and  high-resolution  methods  will  largely  be  ignored.  This  statement  in  effect  reduces 
the  problem  from  that  of  spectral  estimation  to  one  of  parameter  estimation  where 
the  parameters  being  sought  are  horizontal  wavenumbers. 

Using  the  above  statements  as  a  guide,  wavenumber  estimates  can  be  made  by 
identifying  the  peak  locations  of  the  spectral  estimates  for  both  the  2000  m  and  1000 
m  apertures  shown  in  figures  3-2  and  3-3.  The  location  of  the  true  wavenumbers 
are  plotted  as  the  dotted  lines  in  each  of  the  figures.  From  the  plots,  it  is  clear 
the  high-resolution  method  is  capable  of  determining  the  four  discrete  wavenumber 
components  for  both  the  short  and  long  apertures.  In  contrast,  the  classical  result 
is  not  able  to  resolve  the  four  wavenumber  components  for  the  short  aperture  and 
shows  a  bias  in  the  peak  locations  for  the  long  aperture. 

In  the  next  section,  the  high-resolution  method  is  described  and  characterized 
with  emphasis  on  its  ability  to  detect  discrete  wavenumbers  in  the  presence  of  noise. 
The  estimator  will  then  be  used  to  estimate  the  range-dependent  wavenumber  content 
of  the  data  shown  in  figure  3-1. 
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3.2  High- resolution  wavenumber  estimation 


The  problem  of  wavenumber  estimation  for  point-source  acoustic  data  is  intimately 
related  to  the  problem  of  frequency  estimation  from  a  finite  number  of  noisy  discrete¬ 
time  signal  measurements.  For  this  work,  spectral  estimators,  typically  associated 
with  estimating  functions  of  frequency,  are  to  be  used  as  parameter  estimators,  where 
the  parameters  are  horizontal  wavenumbers.  In  this  context,  an  estimator  is  desired 
for  identifying  and  localizing  peaks  in  the  depth-dependent  Green’s  function  that  cor¬ 
respond  to  the  discrete  horizontal  wavenumbers  of  propagating  modes  using  short  data 
apertures.  The  signal  model  is  deterministic  with  added  noise  and  the  requirement 
of  short  data  apertures,  suggests  the  use  of  modern  spectral  estimation  techniques 
for  extraction  of  the  wavenumber  content.  In  particular,  a  method  based  on  au¬ 
toregressive  spectral  analysis  was  chosen  for  its  high-resolution  frequency  estimation 
property  [4]  [71]  [49].  This  method  will  be  investigated  numerically  by  first  applying  it 
to  time-series  data  for  estimating  frequency  and  then  applying  it  to  the  wavenumber 
estimation  problem  of  shallow-water  acoustics.  Where  appropriate,  reference  will  be 
made  to  the  statistical  characteristics  of  the  results  based  on  the  data.  In  particular, 
the  resolution  of  closely  space  frequencies  or  wavenumbers  will  be  investigated  and 
the  variance  of  the  estimates  will  be  considered  for  signals  with  additive  noise.  The 
resolution  will  be  looked  at  with  regard  to  the  theoretical  resolution  of  AR  based 
frequency  estimators  given  by  Marple  [49], 

c.  1-03  .31^ 

^  l)pi’ 

where  5f  is  the  frequency  resolution,  T  is  the  sample  interval  in  seconds,  p  is  AR 
model  order,  and  SNRun  is  the  signal-to-noise  ratio  in  linear  units.  The  variance 
of  AR  spectral  estimators  used  for  frequency  estimation  was  examined  by  Sakai  [71] 
and  found  to  be  proportional  to  the  inverse  of  the  product  of  data  length  and  the 
square  of  SNR,  where  SNR  is  again  given  in  linear  units.  The  results  of  Marple  and 
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Sakai  show  that  both  resolution  and  variance  of  AR  frequency  estimates  depend  on 
noise  level.  Although  the  dependence  on  noise  level  will  not  be  explicitly  studied 
with  respect  to  the  classical  estimators,  it  should  be  clear  that  as  SNR  decreases,  the 
advantage  of  using  the  AR  estimator  will  decrease.  In  particular,  as  SNR  decreases, 
the  denominator  of  (3.1)  decreases  and  the  number  of  AR  model  parameters  has  to  be 
increased  to  achieve  the  same  resolution.  Further,  for  low  enough  SNR,  even  for  high 
model  orders,  the  resolution  will  approach  the  classical  limit  based  on  data  record 
length  and  the  advantages  of  using  an  AR  based  estimator  are  lost.  Therefore,  SNR 
is  a  prime  consideration  when  applying  the  AR  methods  to  real  data. 

This  thesis  work  was  partially  motivated  by  the  work  of  Rajan  and  Bhatta  [65] 
where  other  high-resolution  frequency  estimators  were  studied.  This  work  follows 
their  format  beginning  with  frequency  estimation  for  a  signal  comprised  of  a  single 
tone  plus  noise  and  continuing  on  to  the  acoustic  problem  of  horizontal  wavenumber 
estimation.  The  single- tone  problem  was  treated  by  Rife  and  Boorstyn  [69]  where 
they  derived  the  Cramer-Rao  (CR)  bound  for  the  estimation  of  frequencies  in  a  noisy 
signal.  The  CR  bound  is  the  best  estimate  that  can  be  made  for  a  given  parameter 
using  an  unbiased  estimator  on  the  available  data.  The  signal  model  for  their  analysis 
was  of  the  form, 

k 

=  (3.2) 

i=l 

where  the  sum  is  over  a  finite  number  of  frequencies,  w,  indexed  by  A:,  bt  are  the 
amplitudes,  6i  are  the  phases,  w{t)  is  the  observation  noise,  and  t  is  time.  Assuming 
a  constant  sampling  rate  of  Fg  =  1/T  Hz  with  the  first  sample  taken  at  to,  the 
continuous  signal  model  can  be  written  in  discrete  form  using  the  convention  for 
discrete  time  given  by 

tfi  —  to  fiT  =  (tIo  +  n)T.  (3.3) 
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A  single  sample  is  then  given  at  time  by 

x[n\  =  x{tn)  =  +  w{tn).  (3.4) 

i=l 

For  their  work,  and  the  discussion  that  follows,  the  noise  is  independent  Gaussian 
white  noise  with  zero  mean  and  variance  a^.  When  estimating  a  parameter,  the 
mean-squared  error  is  given  by, 

mse{Q;}  =  S{\oi  —  d|}^  =  var{a}  +  |B(q:)P,  (3.5) 


where  a  is  the  true  parameter  value,  a  is  the  estimate,  S  is  the  expected  value 
operator,  and  B  is  the  bias.  The  CR  bound  for  an  unbiased  estimator,  defined  as  the 
mse  with  B{a)  =  0,  gives  a  minimum  value  for  the  variance  that  can  be  obtained  for 
parameter  estimates  using  the  given  data.  For  a  single-tone.  A:  =  1,  and  known  initial 
phase,  9,  the  CR  bound  is  given  by  [69], 

2 

mse{C)  >  var{a,}  =  427.2, I  2„„p  +  Qy  (3-6) 

where,  b  is  the  amplitude,  rio  is  the  number  of  the  first  sample,  N  is  the  number  of 
samples,  and  P  and  Q  are  defined  as. 


n=0  ^ 


(3.7) 


N{N  -\){2N  -\) 


(3.8) 


n=:0 


This  expression  was  used  as  a  measure  of  comparison  by  Rajan  and  Bhatta  [65]  for 
determining  the  frequency  content  of  noisy  signals  using  eigenanalysis-based  frequency 
estimators.  In  the  next  section,  a  parametric  model  based  AR  spectral  estimator  will 
be  described.  Results  for  estimating  a  single  frequency  using  the  AR  estimator  will 
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be  compared  to  the  CR  bound  above  with  no  =  0  and  with  the  results  of  Rajan 
and  Bhatta.  The  estimator  is  also  applied  to  estimate  temporal  signals  with  two 
closely  spaced  frequency  components  and  the  related  problem  of  estimating  horizontal 
wavenumber  content  for  shallow-water  acoustic  data.  The  multiple  tone  results  are 
compared  to  the  theoretical  results  of  Marple  and  Sakai  along  with  the  numerical 
results  of  Rajan  and  Bhatta  [65]  where  the  above  bound  for  a  single  tone  defined  by 
Rife  and  Boorstyn  does  not  apply. 

3.2.1  Parametric  Modeling  Approach  to  Spectral  Estimation 

Model  based  approaches  to  spectral  estimation  make  use  of  a  parametric  representa¬ 
tion  of  the  data  to  determine  the  spectrum.  It  is  an  approach  based  on  the  assumption 
that  the  measured  data  can  be  related  to  the  input  data  through  a  rational  transfer 
function.  The  data  model  is  represented  by  the  linear  difference  equation, 

q  p 

^kXn-k,  (3.9) 

«=0  Jt=l 

where  {rUn}  is  the  input  driving  sequence  with  associated  parameters  bi,  and  {xn} 
is  the  output  sequence  with  parameters  a*  along  with  their  respective  ^-transforms 
given  by, 

=  Z]  “m-s""",  (3.10) 

m=0 

=  Y.  (3.11) 

m=0 

For  the  given  input  and  output  sequences,  the  transfer  function,  H{z),  is  simply, 

mz)  =  fM.  (3.12) 
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From  these  relationships,  the  power  spectrum  at  the  output  of  a  linear  filter, 
is  related  to  the  input  process,  Pn{z),  by 

p,{z)  =  if(z)H-(iA*)aw  =  (3.13) 

where  *  denotes  the  complex-conjugate.  For  a  white-noise  Gaussian  process  with 
variance  evaluated  along  the  unit  circle,  z  =  exp(z27r/T)  for  — Fs/2  <  /  <  F’s/2, 
Px  becomes 

PmoMf)  =  PAf)  =  (3.14) 

where  A{f)  =  A{exp[i27TfT])  and  B{f)  =  J5(exp  [z27r/r]). 

In  the  given  context,  model  based  spectral  estimation  requires  three  steps  that 
include  selection  of  the  time-series  model,  estimation  of  the  model  parameters,  and 
calculating  the  theoretical  power  spectral  density  (PSD)  from  the  resulting  model. 
The  discussion  above  and  the  derivation  of  model  based  methods  and  other  modern 
estimation  techniques  are  reviewed  by  Kay  and  Marple  [42]  with  practical  applications 
and  computer  algorithms  presented  in  Marple  [49].  In  any  of  these  methods,  the  goal 
is  increased  resolution  and  fidelity  of  frequency  estimates  compared  to  the  classical 
methods.  The  success  of  any  given  method  relies  on  the  ability  to  fit  the  assumed 
model  to  the  data  with  few  parameters.  In  this  work,  a  model  was  sought  that 
could  determine  a  reliable  estimate  of  the  peak  locations  corresponding  to  the  discrete 
wavenumber  spectrum  expected  for  the  depth-dependent  Green’s  function  for  shallow 
water.  This  requirement  is  met  by  a  transfer  function  that  can  be  represented  as  an 
all-pole  model  for  the  spectrum  where  peaks  occur  at  the  pole  locations. 

3.2.2  Autoregressive  Spectral  Estimation 

The  Autoregressive  spectral  estimator  is  defined  by  an  all-pole  model  for  the  rational 
transfer  function  mentioned  above.  For  this  case,  the  first  term  on  the  right  for  the 
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linear-predictor  model  given  by  (3.9)  reduces  to  6o  =  1  and  the  data  model  becomes, 


^  y  ^k^n—k  "b  Wfi. 


(3.15) 


This  model  is  a  strictly  autoregressive  process  of  order  p  where  the  output  sequence 
Xji  is  a  linear  regression  on  itself  with  error  terms  The  present  value  of  the  process 
is  a  weighted  sum  of  past  values  with  a  noise  term.  Using  this  model,  the  spectral 
representation  given  by  (3.14)  reduces  to. 


'ParU)  = 


\A{fW 


(3.16) 


This  equation  can  be  written  explicitly  in  terms  of  the  model  parameters  Ofc  to  show 
the  all-pole  nature  of  the  spectral  estimate,  where  peaks  occur  when  the  denominator 


goes  to  zero. 


Par  = 


(3.17) 


\l  +  T.L,a,exp(-i2TfkTf  '  '  ’ 

Given  the  above  relationship,  an  estimate  of  the  spectrum  is  obtained  by  determining 
the  AR  parameters,  a*  and  This  is  done  formally  by  relating  the  AR  parameters 
to  the  autocorrelation  function  of  the  data  through  the  Yule- Walker  equations  [42]. 
For  the  signal  model,  Xn,  the  autocorrelation  function,  Rxx{k)  is, 

Rxx{k}  ^  ^  '  ^l^n—l+k  "b  Wn+A;^  j  ~  ^  O'lRxxjk 


(3.18) 

The  correlation  between  the  signal  and  noise  in  the  second  term  of  the  last  expression 
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of  (3.18)  can  be  found  by, 


e[Wn+kXn]  =  S 


0,  for  A:  >  0 

for  A:  =  0 


(3.19) 


where  it  was  assumed  that  H{z)  is  a  causal  and  stable  filter  and  Sm  is  the  delta 
function,  with  (5„i  =  1  for  m  =  0  and  0  otherwise.  Imposing  the  asymptotic  limit 
at  z  — >  oo  on  the  transfer  function  gives,  hg  —  lim^^oo^fl-s^)  =  1)  and  yields  the 
Yule- Walker  equations  given  by. 


—  S 


~  Sf=i  o.iRxx{k  1), 
-Ef=l0^^xx(-/)+^T^ 


for  A:  >  0 
for  A:  =  0 


(3.20) 


Given  the  Yule- Walker  equations,  the  AR  parameters  can  be  determined  for  a 
given  model  order  p  by  solving  p  equations  from  (3.20)  for  A:  >  0  and  then  finding 
for  A:  =  0.  In  solving  the  Yule-  Walker  equations  the  model  order  p  is  a  free  parameter 
selected  by  the  user.  Criteria  for  model  order  selection  will  be  discussed  later. 

With  the  model  order  selected,  the  Yule- Walker  equations  can  be  expressed  in 
matrix  form  as 


Rxx{0)  i?ix(-l) 

.Rxi(i)  Rxx{0) 

^ix(p-l)  ■Rxi(p-2) 


/?xx(-(p-l))  ' 

Oi 

’  i?xx(l)  ' 

■Rxi(-(p-2)) 

a2 

—  — 

Rxx{‘^) 

-Rii(o) 

dp 

Rxxip) 

(3.21) 


The  autocorrelation  matrix  above,  R^j;  is  Hermitian  (R^  =  Rix))  Toeplitz  (equal  di¬ 
agonal  elements),  and  is  generally  positive  definite.  Further,  (3.21)  can  be  augmented 
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to  include  the  k  =  0  equation  and  becomes, 


Rxx{0) 

Rxx  (p) 

Using  this  form,  the  AR  parameters  and  a  are  determined  by  solving  (3.22)  with 
p+1  estimated  autocorrelation  estimates,  Rxx{Q), ,  Rxx{p),  and  using  m)  = 

R*xx{fn)- 

An  efficient  method  for  solution  of  (3.22)  is  the  Levinson-Durbin  algorithm  [42]. 
The  algorithm  is  recursive  and  provides  estimates  of  a  number,  p,  AR  parameter  sets 
for  the  selected  model  order  p  and  includes  all  the  lower  model  order  parameter  sets. 
The  final  set  at  order  p  is  the  solution,  with  the  full  complement  of  parameter  sets 
given  by  {an,al},  (021, 0,22,  crj},  {«pi,  ‘ >  Opp,  f^p},  where  the  double  subscript 

denotes  model  order  and  coefficient  number,  respectively.  With  this  notation,  the 
algorithm  is  initialized  by. 


a:xx(-i) 

Rxx{0) 

RxxijP  1) 


Rxxi.  p} 
Rxx{-{p-  1)) 

Rxx{^) 


1 

’(72  ■ 

Oi 

= 

0 

- 1 

_ 1 

1 

0 

1 

(3.22) 


«ii  —  ■~Rxi(i)/Rxx(0) 

af  =  (1  -  |aiiP)Rxx(0) 

and  the  recursion  for  A:  =  2, 3,  •  •  •  is  given  by 

[^•Rxx(^)  T  X/j— 1  Oyk—l,lRxx(,k  l)j  j ^ k—\ 
^ki  l,z  CLkk^k—\^k—i 

al  =  (1-K,|>L,. 


(3.23) 


(3.24) 


Writing  (3.17)  in  terms  of  the  autocorrelation  function  expressed  in  terms  of  the  Yule- 
Walker  equations  and  making  use  of  the  Levinson-Durbin  algorithm  to  extend  the 
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summation  interval  is  the  basis  for  high-resolution  of  the  AR  estimator, 


'Par  = 


a^T 


|1  +  ELi  afcexp(-z27r/A:T)|' 


=  At  ^  exp(— i27r//cT),  (3.25) 


k=:-00 


where 

R:cx[n],  for  |n|  <  p 
^  Efc— 1  ®pfcRix[^  ^])  1^1  ^  P- 

This  compares  to  the  correlogram  method  of  classical  spectral  estimation,  where  the 
summation  is  over  a  finite  number  of  autocorrelation  lags  estimated  from  the  finite 
data  length  [49].  The  effect  of  the  finite  sum  of  the  classical  estimator  is  a  broadening 
of  the  main  peak  of  the  estimate,  as  well  as  the  presence  of  sidelobes  depending 
how  the  data  is  truncated  at  the  endpoints.  This  combination  of  sidelobes  and  peak 
broadening  for  the  classical  estimator  limit  the  frequency  resolution  for  short  aperture 
data.  The  AR  spectral  estimator  is  not  limited  by  windowing  effects  on  the  data  and 
resolution  is  dependent  on  model  order  selection  and  data  aperture. 

The  solution  of  the  AR  model  parameters  can  be  recast  in  terms  of  the  theory  of 
linear  prediction.  Inspection  of  equation  (3.18)  shows  it  to  be  the  equation  of  linear 
prediction  where  x  is  to  be  predicted  given  p  previous  samples  of  the  signal  [42]. 


Rxx 


x  =  -J2akXn-k,  (3.27) 

In  this  case,  the  best  linear  prediction  for  x  will  minimize  the  prediction  error  power, 
Qp,  where 

Qp  =  S[\Xn  -  Xnf]-  (3.28) 

It  can  be  shown  that  0;*,  =  CLpk  for  ^  =  T" '  ^^od  Qpm.n  —  ^p  ^^e  best  linear 
predictor  for  x  is, 

p 

X  —  ^  ^  RpkXjx-k-  (3.29) 

k=\ 

The  above  states  that  the  parameters  {0^1,0^27  •  •  •  ,o,kk},  (where,  as  before,  the  sub- 
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scripts  k  and  p  refer  to  the  order  and  index)  along  with  represent  the  kth  —  order 
linear  predictor  and  minimum  prediction  error  power  and  are  identical  to  AR  pa¬ 
rameter  estimation  allowing  the  theory  for  one  to  apply  to  the  other.  Thus,  linear 
prediction  theory  can  be  used  for  estimating  AR  parameters.  This  is  very  powerful  in 
that  it  allows  both  forward  and  backward  prediction  schemes  to  be  used  in  estimating 
AR  parameters.  Writing  the  pth  —  order  forward  prediction  error  cTp„  as  ep„  given  by, 


+ 

2^k—\  ^pk^n—k 

Xn 

+ 

^k=li^P~^,k  ^pp^p—l^p—k)^Ti—k  “1"  (^pp^n~p 

(3.30) 

+ 

^ppbp-l^n-l) 

where 

p 

6pn  =  ^n—p  T  ^  !  ^pk^n—p+kj  (3.31) 

k=\ 

is  the  backward  prediction  error,  and  the  term  is  predicted  from  the  samples 
x„_p+i,  •  •  • ,  Xn,  it  can  be  shown  that  the  backward  prediction  coefficients,  6p„,  are 
complex  conjugates  of  the  forward  prediction  coefficients,  ap„  [49]. 

Having  established  the  relationship  between  AR  parameters  and  the  forward  and 
backward  prediction  coefficients,  a  method  for  estimation  was  chosen  that  combines 
the  forward  and  backward  linear  prediction  algorithms.  By  combining  the  error  statis¬ 
tics  of  both  the  forward  and  backward  predictors,  more  error  points  are  generated 
resulting  in  an  improved  estimate  of  the  AR  parameters.  The  method  of  solution  is 
that  of  Marple  [49]  and  is  a  forward  and  backward  least  squares  approach  referred  to 
as  the  modified  covariance  method. 

Writing  the  backward  prediction  error,  bpn,  of  equation  (3.31)  as  the  vector,  e*. 
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and  the  forward  error,  ep„,  as,  e^,  yields  the  two  {N  —  p)  length  vectors  given  by, 


(3.32) 


where  the  double  subscript  has  been  replaced  by  a  single  subscript  for  the  order  and 
sample  number  indices  are  given  explicitly.  The  combined  error  can  be  expressed  as 
the  length  2{N  —  1)  vector.  Bp, 


(3.33) 


(3.34) 


(3.35) 


70 


The  linear  prediction  coefficient  vector  is  given  by, 


a^**  = 

“p 


(3.36) 


Using  the  above,  a  least  squares  approach  is  sought  which  minimizes  the  average 
of  the  forward  and  backward  linear  prediction  squared  errors  over  the  available  data, 

Pp  =  \  E  E  =^efep  =  ^[(e^)"e^  +  (ej)^ej].  (3.37) 


n=pH-l 


n— p+1 


This  leads  to  a  set  of  normal  equations  to  be  solved  as  given  by 


Rpa;' 


(3.38) 


with. 


t"t,  +  jtJt;j, 


(3.39) 


and  Op  a  p-element  zero  vector.  The  individual  elements  of  Rp  are  given  by. 


^p[bi]  =  E  {x*[n-l]x[n- j]-\-x[n--p-\-l]x*[n-p-\r  j]) ,  for  0  <  jf ,  A:  <  p. 

n=p+l 

(3.40) 

Marple  [49]  describes  an  efficient  algorithm  for  the  solution  of  (3.38)  with  a  couple 
of  conditions.  One  requirement  is  that  for  Rp  to  be  non-singular,  p  <  2N/3,  i.e.,  the 
order  selected  can  be  no  larger  then  2/3  the  data  length.  Another  caveat  is  that  for  a 
purely  harmonic  signal  the  recursion  of  the  Levinson-Durbin  algorithm  will  terminate 
for  any  k  where  \akk\  =  1  since  =  0.  Thus,  for  signals  comprised  of  pure  or  near 
sinusoids,  the  recursion  algorithm  may  terminate  before  reaching  the  specified  model 
order. 
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The  modified-covariance  AR  algorithm  will  be  implemented  on  a  signal  with 
noise  for  estimating  discrete  frequency  components  and  compared  with  the  high- 
resolution  estimation  results  of  Rajan  and  Bhatta  [65]  and  the  CR  bound  of  Rife 
and  Boorstyn  [69].  It  will  then  be  used  to  estimate  wavenumbers  of  the  propagating 
modes  for  shallow  water  acoustic  fields. 


3.3  Performance  of  Modified  Covariance  AR  fre¬ 
quency  estimator 

In  order  to  judge  the  performance  characteristic  of  the  AR  algorithm  to  estimate 
wavenumbers,  several  tests  were  performed  on  time  series  data  for  a  signal  plus  noise. 
These  tests  were  adopted  from  the  paper  of  Rajan  and  Bhatta  [65]  for  comparison 
with  their  results.  The  time  domain  signals  were  chosen  with  frequency  and  time¬ 
sampling  intervals  that  are  typical  of  wavenumber  differences  and  range-sampling 
intervals  for  acoustic  measurements  made  in  shallow-water  waveguides.  The  signal 
data  for  the  numerical  studies  were  generated  using 

L 

yW  =  A;  exp[f27r /^nT]  +  w[n],  n  —  —  (3-41) 

i=l 

where  Aj  are  the  amplitudes,  fi  are  the  frequencies,  T  is  the  sampling  interval  in 
seconds,  and  w[n]  is  the  added  noise  assumed  to  be  zero-mean  white  Gaussian.  For 
all  studies,  the  amplitudes  were  assumed  equal  with  a  value  of  1,  and  the  sampling 
interval  was  1  second. 

3.3.1  Frequency  Estimation  in  the  Presence  of  Noise 

The  AR  algorithm  was  used  to  estimate  the  frequency  content  of  complex  sinusoids 
with  added  noise.  Frequency  was  determined  by  the  peak  location  of  the  estimated 
signal  spectrum.  For  the  single  frequency  case,  the  sum  in  (3.41)  reduces  to  a  single 
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SNR  (dB) 

MUSIC 

ESPRIT 

AR 

C-R 

5 

-51.06 

-51.78 

-  58.14 

-58.16 

10 

-56.56 

-58.10 

-63.66 

-63.16 

20 

-66.52 

-68.50 

-73.94 

-73.16 

30 

-76.38 

-78.54 

-82.82 

-83.16 

40 

-86.20 

-88.56 

-90.86 

-93.16 

Table  3.1:  SNR  (dB)  vs  lOlogio(n^se)  for  estimation  of  frequency  =  2.3/27r  Hz. 

term  and  the  spectrum  will  contain  a  single  peak.  Estimates  were  made  for  increasing 
noise  levels  as  indicated  by  the  signal  to  noise  ratio  (SNR).  SNR  for  this  study  is 
defined  as  10/o^io(crf/cr^),  where  and  cr^  are  the  mean-square  signal  power  and 
noise  power,  respectively.  Estimates  were  made  for  a  fixed  data  length  of  40  samples 
(N  =  40),  and  100  realizations  of  the  additive  noise  were  generated  to  determine 
performance  statistics.  AR  model  order  was  determined  by  trial  and  error  at  each 
SNR  as  the  maximum  order  that  gave  only  a  single  peak  for  all  realizations  of  the 
signal.  For  SNR  <  30  dB  a  model  order  of  13  was  used  which  is  approximately  equal 
to  N/3.  The  mean  squared  error  of  (3.5)  was  used  as  a  basis  for  performance  and 
compared  to  the  CR  bound  given  by  (3.6).  Results  are  given  in  dB  calculated  as 
lOlogio(mse).  The  ideal  result  for  an  unbiased  estimator  is  one  that  achieves  the  CR 
bound.  In  addition  to  the  CR  bound,  the  AR  results  are  compared  with  the  estimates 
of  the  MUSIC  (Multiple  Signal  Classification  Method)  and  ESPRIT  (Estimation  of 
Signal  Parameters  via  Rotational  Invariance  Techniques)  determined  by  Rajan  [65], 
where  his  results  have  been  converted  to  match  the  scale  defined  for  this  study. 

Frequency  estimates  were  made  for  individual  frequencies  /i  =  2.5/27r  Hz  and 
/2  =  2.3/27r  Hz  for  SNR  values  from  5  to  40  dB.  The  mse  results  for  estimating  /2 
are  shown  in  table  3.1.  Bias  was  calculated  for  each  of  the  two  frequencies  at  the 
different  SNR  and  are  shown  in  table  3.2.  From  the  results  given  in  tables  1  and  2 
the  AR  estimator  shows  the  following  behavior. 

1.)  The  mean-squared  error  is  better  for  the  AR  estimator  then  either  MUSIC  or 
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SNR  (dB) 

freq  =  2.5/  tt  Hz 

freq  =  2.3  /  tt  Hz 

5 

-2.56x10“'* 

10 

7.88x10“* 

20 

-7.38x10-* 

-1.02.xl0“* 

30 

-4.06x10“* 

40 

-2.39x10“* 

-3.47x10-5 

Table  3.2:  SNR  (dB)  vs  bias  of  estimates  for  frequencies  fi  =  Hz  and  /a  = 

2.3/27r  Hz. 

ESPRIT  and  approaches  the  CR  bound  at  all  SNR  levels. 

2.)  The  bias  error  is  of  the  same  order  or  better  then  MUSIC  and  ESPRIT  at  all 
SNR  levels. 

Although  these  results  are  compared  with  the  CR  bound,  the  estimation  results 
at  each  different  SNR  had  a  small  bias  which  violates  the  criterion  of  an  unbiased 
estimator.  Given  a  limited  number  of  realizations,  this  discrepancy  can  cause  the 
estimator  to  perform  better  than  the  bound  as  the  results  in  table  3.1  suggest  at 
low  SNR.  Nevertheless,  the  comparision  is  still  of  use  for  evaluating  results  from 
different  estimators.  Further,  the  bias  for  SNR  i  10  dB  was  sufficiently  small  that  mse 
effectively  gives  a  measure  of  the  variance  of  the  estimates.  The  variance  of  5.3a;10''^ 

I 

at  30  dB  SNR  also  compares  favorably  to  the  estimate  2.5x10“*  as  determined  by 
Sakai’s  [71]  result  which  is  proportional  to  SNRj^/N .  Finally,  it  should  be  noted, 
in  order  to  measure  the  small  variances  in  the  estimated  frequencies  at  high  SNR,  a 
large  number  of  terms  was  used  in  the  summation  interval  of  (3.25)  to  give  sufficient 
discritization  in  the  frequency  domain.  For  SNR  of  30  dB,  2^*  points  were  used  in 
the  summation. 

3.3.2  Frequency  Resolution 

The  next  study  uses  the  AR  estimator  to  identify  two  closely  spaced  frequencies 
in  the  presence  of  added  noise.  This  is  representative  of  the  case  of  closely  spaced 
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SNR  (dB) 

ESPRIT 

AR 

20 

-53.74 

-61.16 

30 

-63.66 

-70.28 

40 

-73.74 

-77.70 

50 

-83.72 

-79.48 

Table  3.3:  SNR  (dB)  vs  lOlogio(mse)  for  estimation  of  fo  using  ESPRIT  and  AR. 
fo  =  0A  Hz,  5/ =  0.011  Hz. 

wavenumbers  for  a  shallow  water  waveguide.  The  signal  model  of  (3.41)  was  used 
with  two  frequencies  of  various  separation  frequencies,  5 /,  that  ranged  from  0.009  Hz 
to  0.027  Hz.  Using  equation  (3.1),  at  20  dB  SNR,  a  minimum  model  order  of  13  would 
be  required  to  resolve  5f  =  0.009  Hz.  The  frequencies  used  to  generate  the  complex 
signal  were  fi  =  fo  =  0.4Hz  and  /2  =  fo  ~  Hz.  Each  signal  was  comprised  of  40 
samples  and  100  realizations  were  made  of  the  noise  at  each  SNR.  Model  order  was 
again  selected  by  trial  and  error  to  yield  two  spectral  peaks  for  all  signal  realizations 
with  a  fixed  order  of  13  for  SNR  under  40  dB.  Table  3.3  gives  the  mse  estimates 
vs  SNR  for  estimation  of  /„  =  0.4  Hz  and  a  frequency  separation  of  0.011  Hz.  In 
these  cases  where  both  estimators  were  able  to  resolve  the  two  frequencies,  the  results 
indicate  that  the  AR  estimator  performs  better  than  ESPRIT  at  all  SNR  for  resolving 
closely  spaced  frequencies.  Similar  results  were  observed  for  the  other  estimates  of 
closely  spaced  frequencies. 

3.3.3  Size  of  Data  Vector  and  Model  Order  Selection 

In  the  previous  studies,  AR  model  order  was  selected  by  trial  and  error  based  on 
the  number  of  known  sinusoids  in  the  signal.  However,  the  resolution  and  variance 
of  the  estimator  were  shown  by  Marple  [49]  and  Sakai  [71]  to  be  a  function  of  total 
data  length,  model  order,  and  SNR.  Further,  although  the  all  pole  model  of  the  AR 
estimator  suggest  that  model  order  be  selected  equal  to  the  number  of  sinusoids  in 
the  signal,  the  presence  of  noise  requires  a  higher  model  order  to  best  fit  the  model  to 
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Samples 

ESPRIT 
%  correct 

AR 

%  correct 

40 

20 

34 

60 

21 

55 

80 

88 

81 

100 

99 

98 

Table  3.4:  Probability  of  success  for  estimation  of  f  =  0.4  Hz  with  SNR  =  10  dB, 
model  order  p  =  20,  and  increasing  data  length.  Percent  correct  for  estimating  (fl  ± 
0.001)  Hz. 

the  data.  As  noise  increases,  the  number  of  model  parameters  required  to  best  fit  the 
data  also  increases.  In  this  section,  the  effect  of  data  length  is  examined  as  it  effects 
resolution  in  resolving  signals  comprised  of  closely  spaced  sinusoids.  An  additional 
study  is  then  run  to  examine  the  effect  of  model  order  selection  on  resolution  and 
variance  and  the  presence  of  spurious  peaks  in  the  spectral  estimates. 

The  first  of  these  tests  was  conducted  to  identify  two  frequencies,  /i  =  0.4  Hz  and 
/2  =  0.39  Hz,  with  SNR  =  10  dB  for  various  data  lengths.  Data  length  was  increased 
from  40  points  to  100  points,  and  the  probability  of  success  in  identifying  /i  ±  0.001 
Hz  was  determined  for  100  realizations  of  the  noisy  signal.  For  SNR  =  10  dB  and 
N  =  40,  a  model  order  of  20  was  necessary  to  theoretically  resolve  the  two  freqencies. 
This  was  used  for  all  data  lengths  in  the  study.  The  results  are  shown  in  table  3.4  and 
compare  reasonably  well  with  the  ESPRIT  results  of  Rajan  and  Bhatta  [65].  Overall, 
the  AR  estimator  performs  better  then  MUSIC  and  ESPRIT  for  short  data  lengths 
of  40  and  60,  and  on  par  with  them  for  longer  data  lengths. 

A  similar  study  was  conducted  to  examine  the  effect  of  model  order  selection  on 
frequency  estimates.  The  signal  was  comprised  of  40  samples  with  two  sinusiods  at 
fl  =  OAHz  and  /2  =  0.3bHz  with  an  SNR  of  20  dB.  Of  particular  interest  is  the  effect 
of  resolution  on  model  order,  and  also  the  presence  of  spurious  peaks  in  the  spectrum 
for  high  model  orders.  The  effect  of  model  order  selection  on  the  spectral  estimates  is 
clearly  shown  in  figures  3-4  and  3-5.  In  each  of  the  figures,  the  left  column  shows  the 
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1 


ARSE(4) 


ARSE(4) 


Frequency  (Hz)  Frequency  (Hz) 

Figure  3-4:  Spectral  estimates  for  100  realizations  of  signal  with  2  sinusoids  in  noise. 
SNR  =  20  dB,  fi  =  QAHz,  f2  =  O.SbHz.  AR  model  order  is  indicated  at  top  of  each 
plot. 


results  of  each  spectrum  for  the  100  signal  realizations  plotted  on  top  on  one  another. 
The  right  column  shows  the  spectra  for  all  realizations  lined  up  next  to  one  another 
with  any  peaks  greater  than  60  dB  down  indicated  in  white.  This  display  gives  an 
indication  of  both  peak  width  and  location.  It  also  illustrates  the  stability  of  the  peak 
locations  for  the  two  excited  frequencies  for  all  realizations.  This  type  of  behavior 
is  typical  of  wavenumber  estimates  for  range-independent  shallow  water  waveguides 
and  is  used  to  identify  the  prominant  propagating  modes  using  the  sliding  window 
transform  method.  From  the  plots,  it  is  clear  that  as  model  order  increases,  the 
peak  widths  decrease.  However,  as  model  order  increases  past  N/3,  spurious  peaks 
can  be  introduced  into  the  spectra.  The  presence  of  extra  peaks  causes  ambiguity 
in  the  selection  of  the  peaks  which  represent  the  true  signal.  Model  order  selection 
for  AR  spectral  estimates  with  no  a  priori  knowledge  must  strike  a  balance  between 
resolution  and  the  introduction  of  spurious  peaks. 

Marple  [49]  reviews  several  methods  and  criteria  put  forth  for  determining  model 
order  selection,  including  the  final  predicion  error  (FPE)  [1]  and  Akaike  information 
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Figure  3-5:  Spectral  estimates  for  100  realizations  of  signal  with  2  sinusoids  in  noise. 
SNR  =  20  dB,  /i  =  OAHz,  /2  =  0.35Hz.  AR  model  order  is  indicated  at  top  of  each 
plot. 


criterion  (AIC)  of  Akaike  [2],  where  Gaussian  noise  statistics  are  assumed.  Each  of 
these  criteria  seek  a  minimum  value  in  the  prediction  error  variance  with  a  penalty 
applied  for  increased  model  order  p  after  the  minimum  is  reached.  Marple  also  de¬ 
scribes  the  minimum  description  length  (MDL)  and  criterion  autoregressive  transfer 
(CAT)  function.  After  examining  all  these  criteria,  Marple  states  that  the  results 
of  applying  each  to  real  data  are  inconsistent.  Various  empirical  studies  are  then 
presented  which  suggests  that  for  noisy  data  a  model  order  of  N/3  to  N/2  yield  good 
results  on  short  data  segments.  This  hypothesis  is  tested  using  the  AR  estimator  used 
for  this  work.  For  this  study  40  data  points  were  used  with  an  SNR  of  20  dB  for  the 
same  frequencies  used  in  the  previous  case,  /i  =  0.4  Hz  and  /2  =  0.39  Hz.  Success 
rates  were  calculated  for  100  realizations  of  the  noise  for  model  order  of  between  2 
and  25,  with  model  order  p  =  13  approximately  N/3.  The  results  shown  in  table  3.5 
confirm  the  empirical  results  cited  by  Marple  as  well  as  those  shown  in  figures  3-4  and 
3-5.  Success  rates  are  maximized  for  a  model  order  p  =  13.  Results  marked  by  an 
X,  indicate  that  no  peaks  were  observed  in  that  bin  for  any  realizations.  It  was  ob- 
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Samples 

mse(fl) 

%  (fl  ±  0.001) 

mse(f2) 

%  (f2  ±  0.001 

2 

X 

X 

0 

1 

5 

X 

X 

0 

0 

10 

-57.92 

65 

-58.98 

71 

13 

-62.64 

78 

-64.14 

88 

15 

-62.48 

83 

-63.42 

83 

20* 

-22/-64 

79 

-26/-66 

89 

25 

-18 

59 

-28 

79 

Table  3.5:  Probability  of  success  for  estimation  of  f  =  0.4  Hz  with  SNR  =  10  dB  and 
increasing  model  order. 

served  that  as  model  order  was  increased  beyond  N/Z,  spurious  peaks  contaminated 
the  spectra.  However,  for  the  case  of  p  =  20,  where  two  values  are  listed,  a  small 
number  of  outliers  were  removed  using  a  priori  information  and  success  rates  were 
improved.  However,  outlier  removal  requires  information  that  may  not  be  available 
so  is  not  generally  a  viable  approach  to  improving  on  estimates. 

A  final  comment  about  model  order  selection  is  in  regard  to  the  true  nature  of 
the  process  being  modelled  by  the  AR  parameters.  For  a  true  AR  model,  there  exists 
a  given  model  order  p  for  which  the  prediction  error  variance  estimate  is  minimized. 
For  model  orders  greater  then  p,  the  prediction  error  variance  is  flat  and  and  gives 
a  method  for  selecting  the  exact  order  of  the  AR  process.  However,  for  a  non-AR 
process,  prediction  error  variance  is  a  monontonically  decreasing  function  of  p  and  no 
minimum  will  be  reached.  This  contributes  to  the  problem  of  selecting  a  model  order 
a  priori  for  a  non-AR  process. 

3.3.4  Application  of  AR  estimator  to  acoustic  data 

It  is  now  of  interest  to  use  the  AR  spectral  estimator  to  determine  the  horizontal 
wavenumbers  for  the  normal  modes  propagating  in  a  shallow- water  waveguide.  Cor¬ 
recting  for  geometric  spreading  and  adding  a  noise  term,  the  signal  to  be  analyzed 
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Depth  (m) 

c,  (m/s) 

P  (g  /cm^) 

a  (nepers/m) 

0-40 

1515.0 

1.0 

0.0 

40+ 

1800.0 

1.56 

0.001 

Table  3.6:  Pekeris  ocean  model.  Frequency  =  100  Hz,  source  depth  =  20.5  m  and 
receiver  depth  =  15.5  m. 


Mode  No. 

kn  (1/m) 

1 

0.409260 

0.179 

0.204 

2 

0.392023 

0.181 

0.074 

3 

0.361364 

-0.003 

-0.178 

Table  3.7:  Modal  eigenvalues  and  corresponding  amplitudes  at  source  and  receiver 
for  Pekeris  ocean  model. 

can  be  represented  in  terms  of  the  acoustic  pressure,  P(nAr),  by  [65] 

L 

y[n]  =  VnArP[nAr]  +  rz;[nAr]  =  ^  Ai'^i  +  w[nAr],  (3.42) 

/=i 

where  the  range  interval  is  sampled  by  Ar,  and  is  the  mode  function  with 

horizontal  wavenumber  ki  and  amplitude  Ai.  The  objective  is  to  estimate  the  ki  from 
the  data  measured  along  r.  To  generate  data,  KRAKEN  [59],  a  normal-mode  acoustic 
propagation  code,  was  used.  The  output  of  KRAKEN  gives  complex  pressure  fields 
as  a  function  of  range  as  well  as  eigenvalue  information  for  the  input  environment 
and  source  frequency.  The  environment  model  used  was  a  Pekeris  ocean  model  with 
parameters  listed  in  table  3.6.  The  source  frequency  was  100  Hz,  with  source  and 
receiver  at  depths  of  20.5  m  and  15.5  m,  respectively.  The  field  was  sampled  on  a  1  m 
range  grid.  Mode  numbers  and  modal  amplitudes  at  the  source  and  the  receiver  are 
listed  in  table  3.7.  The  KRAKEN  results  show  that  mode  3  is  only  weakly  excited 
at  the  source  depth  and  will  be  difficult  to  estimate.  For  the  given  field,  wavenumber 
estimates  are  made  and  mse  calculated  for  various  aperture  lengths  and  SNR  levels. 
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SNR  (dB) 

Aperture  (m) 

Mode  No. 

ESPRIT 

AR 

40 

240 

1 

-81.32 

-83.00 

2 

-69.74 

-72.46 

3 

-19.38 

-39.68 

40 

360 

1 

-88.61 

-83.00 

2 

-77.58 

-78.70 

3 

-53.84 

-61.62 

20 

360 

1 

-69.40 

-73.98 

2 

-61.76 

-66.82 

20 

480 

1 

-81.14 

-79.50 

2 

-71.62 

-71.48 

Table  3.8:  SNR  (dB)  vs  lOlogio(mse)  of  eigenvalues  for  different  data  aperture 
lengths. 

AR  model  order  was  selected  as  N/3  as  determined  by  total  aperture  length.  Table 
3.8  shows  the  results  of  these  calculations.  The  results,  as  indicated  by  small  mse, 
show  that  the  strong  modes  are  detected  for  all  aperture  lengths  and  SNR  levels.  The 
weak  third  mode  is  identified  with  confidence  only  for  high  SNR  and  long  aperture. 
Consistent  with  the  previous  study  of  time  series  data,  the  AR  estimator  performs  as 
well  as  or  slightly  better  than  the  other  high-resolution  methods. 

3.4  Wavenumber  Estimates  for  Range-Dependent 
Synthetic  Acoustic  Data 

Rajan  and  Bhatta  [65]  also  studied  wavenumber  estimation  using  high-resolution  esti¬ 
mators  for  range-dependent  acoustic  environments.  Their  method  examined  wavenum¬ 
ber  estimates  resulting  from  short  aperture  transforms  using  data  from  a  model  that 
included  bathymetry  with  a  constant  slope  in  range.  The  bathymetric  slope  was  1.67 
m/km  ,  which  equates  to  an  angle  of  0.955°.  In  general,  their  wavenumber  estima¬ 
tors  yielded  wavenumber  values  that  were  between  those  expected  for  the  deepest 
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and  shallowest  portion  of  the  aperture.  They  concluded  that  errors  associated  with 
wavenumber  estimates  under  the  assumption  of  locally  range-independent  segments 
are  small.  Further,  given  the  previous  studies,  the  errors  associated  with  this  slight 
range  dependence  would  be  indistinguishable  in  the  presence  of  noise.  The  behavior 
described  by  Rajan  and  Bhatta  was  observed  for  wavenumber  estimates  made  on 
an  independent  data  set  provided  as  part  of  a  recent  inversion  techniques  workshop 
(ITW)  [40]. 

3.4.1  Inversion  Techniques  Workshop 

The  ITW  was  conducted  in  order  to  evaluate  the  current  state  of  the  art  in  geoacous¬ 
tic  inversion  algorithms  for  application  in  range-dependent  shallow-water  waveguides. 
Synthetic  acoustic  pressure  field  data  were  generated  for  several  test  case  environ¬ 
ments  and  given  to  participants  for  analysis.  The  data  were  provided  at  a  number  of 
discrete  frequencies  on  both  horizontal  and  vertical  arrays.  Properties  of  the  water 
column  were  given,  as  were  source  and  receiver  depths.  Bathymetry  was  assumed 
known  to  within  a  meter.  No  information  was  given  regarding  the  sediment  proper¬ 
ties,  and  participants  were  asked  to  invert  for  compressional  wave  speed,  attenuation, 
and  density  in  the  bottom.  A  total  of  four  different  synthetic  tests  cases  (TC0-TC3) 
were  given  for  analysis.  TCO  was  a  calibration  case,  shown  in  figure  3-6,  and  was 
characterized  by  a  constant  bathymetric  slope  where  the  stratigraphy  in  the  bottom 
followed  the  bathymetry.  The  input  model  for  TCO  was  given  so  that  participants 
could  calibrate  their  forward  propagation  models  to  the  given  synthetic  data.  TCO 
was  used  in  this  work  to  examine  the  characteristic  of  the  AR  wavenumber  estimator 
for  a  range-dependent  problem.  Because  the  model  was  given,  KRAKEN  could  be 
used  to  calculate  the  range-dependence  of  the  eigenvalues.  Sliding  window  AR  esti¬ 
mates  could  be  made  and  wavenumber  estimates  compared  to  the  expected  values  at 
a  given  range  step.  The  sliding  window  estimator  results  are  plotted  in  figure  3-7  for 
a  1000  m  aperture  with  a  step  length  of  50  m.  The  data  was  provided  on  a  5  meter 
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TCO  -  Calibration 


Figure  3-6:  Environment  model  for  Inversion  Workshop  TCO. 


range  grid  resulting  in  an  N/3  model  order  selection  of  66  for  a  1000  m  aperture.  The 
estimates  for  the  50  Hz  data  and  the  true  values  are  shown  in  figure  3-8.  The  true 
values  are  bounded  by  the  values  determined  for  the  endpoints  of  the  1000  m  aperture 
and  represent  the  minimum  and  maximum  values  within  the  aperture.  Similar  to  the 
Rajan  study,  the  AR  estimated  wavenumbers  fall  between  the  minimum  and  maxi¬ 
mum  values  for  most  range  intervals.  And,  although  the  resulting  wavenumber  values 
did  not  depart  very  much  from  expected  values,  the  inconsistent  prediction  behavior, 
shown  as  oscillations  about  the  true  values,  is  undesirable.  One  can  conclude  from  this 
study  that  the  locally  range-independent  assumption  cannot  be  applied  for  constant 
slopes  in  bathymetry  of  even  1  degree,  and  that  bathymetry  must  be  accounted  for 
in  the  wavenumber  estimates.  Harrison  and  Siderius  [37]  suggest  a  method  for  using 
range-independent  models  to  calculate  range-dependent  effects  through  the  use  of  an 
effective  geometry.  In  their  work,  an  environment  similar  to  TCO  is  considered  and 
some  of  their  ideas  may  be  applicable  to  determining  wavenumbers  for  the  continu¬ 
ously  varying  waveguide.  Their  ideas  are  similar  in  nature  to  those  of  Fernandez  [24], 
where  horizontal  wavenumber  analysis  is  investigated  for  range-dependent  environ- 
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tcO  (50  Hz)  1000  X  50  m  AR 


Figure  3-7:  Sliding  Window  spectral  estimates  for  50  Hz  TCO.  Triangles  are  the 
range-dependent  horizontal  wavenumbers  determined  using  KRAKEN  plotted  over 
the  wavenumber  estimates. 


TcO  50  Hz 


Figure  3-8:  Wavenumber  estimates  for  TCO  (dots)  plotted  with  values  for  the  starting 
range,  center  range,  and  end  range  of  the  transform  aperture. 
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TC2  -  Shelfbreak 


◄ -  5000m  - ► 

Figure  3-9:  Shelfbreak  environmental  model  for  Inversion  Workshop  TC2. 


ments  by  migrating  the  measured  field  to  an  effective  range  grid  to  approximate  a 
range-independent  environment.  Another  area  of  investigation  for  wavenumber  anal¬ 
ysis  in  range-dependent  environments  would  be  the  application  of  time-varying  au¬ 
toregressive  estimators  (TVAR)  [43]  [82].  The  TVAR  model  allows  for  time- variation 
in  the  AR  model  coefficients  for  non-stationary  time  series  data  that  can  include  both 
trends  and  discrete  changes  in  the  signal.  However,  as  the  AR  estimator  discussed 
here  does  not  perform  well  for  the  continuously  varying  environment,  TCO  and  TCI 
(both  constant  slope)  were  not  analyzed  for  the  workshop.  Instead,  emphasis  was 
placed  on  examining  wavenumber  estimates  for  range-dependent  cases  where  bathy¬ 
metric  effects  are  minimized,  this  includes  regions  of  slightly  undulating  bathymetry, 
and  regions  of  constant  local  bathymetry.  An  example  of  the  later  is  TC2,  shown  in 
figure  3-9,  where  a  fiat  bottom  occurs  after  a  shelf-break.  The  wavenumber  estimates 
at  50  Hz  for  TC2  are  shown  in  figure  3-10  and  appear  to  be  range-independent  for 
ranges  after  the  break.  Workshop  organizers  confirmed  this  to  be  true.  In  this  case, 
the  data  again  was  sampled  on  a  5  meter  grid,  and  a  1000  m  sliding  window  estimator 
was  used  with  a  step  size  of  50  m  and  a  model  order  of  20.  The  low  model  order 
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was  used  because  the  data  is  noise  free  and  approximates  a  sum  of  pure  sinusoids 
after  the  shelf  break  causing  the  Levinson-Durbin  algorithm  to  terminate  at  a  low 
model  order.  The  effect  of  the  low'  model  order  in  the  range-dependent  region  is  low 
resolution  as  shown  in  the  figure.  However,  if  model  order  is  selected  high  enough 
in  the  range-dependent  region,  results  similar  to  those  of  TCO  could  be  expected. 
By  using  an  algorithm  with  an  adaptive  model  order  selector  as  the  window  stepped 
through  range,  better  resolution  of  the  modes  in  the  range-dependent  region  could  be 
obtained  as  shown  in  figure  3-11.  For  this  result,  model  order  was  linearly  decreased 
from  p  =  66  to  p  =  20  as  the  window  began  to  overlap  the  range-independent  region. 
This  result  highlights  the  effect  of  model  order  selection  on  the  AR  estimates  as  dis¬ 
cussed  previously.  Additionally,  sliding  window  wavenumber  estimates  were  applied 
to  the  TC2  data  with  40  dB  additive  white  noise  added.  In  this  case,  model  order 
was  fixed  at  66  and  a  1000  m  sliding  window  aperture  used.  The  results  are  similar 
to  the  results  for  the  adaptive  estimator  with  no  noise.  The  effect  of  the  noise  on  the 
estimates  were  minimal,  but  allowed  the  recursion  in  the  AR  algorithm  to  finish  at 
all  ranges  eliminating  the  requirement  for  adaptive  model  order  selection. 
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Figure  3-10:  Wavenumber  estimates  showing  range-independent  region  after  the  shelf- 
break  for  Inversion  Workshop  TC2. 


TC2  k^(r)  (50  Hz)  1000  m  AR 


Figure  3-11:  Wavenumber  estimates  using  adaptive  AR  model  order  selector  showing 
wavenumber  evolution  for  Inversion  Workshop  TC2. 
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TC2  k^(r)  (50  Hz)  1000  m  AR 


Figure  3-12:  Wavenumber  estimates  for  order  66  AR  estimator  showing  wavenumber 
evolution  Inversion  Workshop  TC2  data  with  60  dB  noise  added. 
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TC3  -  Flat  Bottom 


◄ -  5000m  - ► 


Figure  3-13:  Environment  model  for  Inversion  Workshop  TC3. 


3.4.2  Range-Dependent  Sediment  -  TC3 

A  particularly  interesting  data  set  for  analysis  was  TC3  where  the  bathymetry  was 
flat  and  range-dependence  occurred  solely  in  the  sediment.  The  geometry  and  other 
known  parameters  are  shown  in  figure  3-13,  along  with  50  Hz  data  on  a  5  meter 
range  grid  for  the  upper  receiver  in  figure  3-1.  No  noise  was  added  to  the  data.  The 
only  a  priori  information  given  about  the  bottom  is  that  “there  is  an  intrusion  in 
the  sediment”  [40].  No  other  information  about  the  intrusion  is  given.  This  point 
is  important  because  depending  on  the  inversion  algorithm  used,  assumptions  may 
have  to  be  made  about  the  nature  of  the  intrusion.  Features  such  as  location  and 
shape  may  have  to  be  assumed  [13]  that  have  an  impact  on  the  inversion  results. 
Techniques  based  on  wavenumber  estimates  need  no  such  assumption  and  are  only 
limited  by  their  ability  to  resolve  features  using  short  data  apertures. 

Using  the  sliding  window  AR  estimator  on  the  50  Hz  data  with  a  model  order 
of  30  yields  the  wavenumber- range  plot  shown  in  figure  3-14.  From  the  figure  two 
range-independent  regions  can  be  identified  directly  by  the  constant  wavenumber 
values.  The  other  regions  of  low  resolution  are  a  result  of  the  fixed  model  order 
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TC3  k^(r)  (50  Hz)  1000  m  AR 


Figure  3-14:  Wavenumber  estimates. for  range-dependent  sediment  of  Inversion  Work¬ 
shop  TC3. 

used  with  the  sliding  window  estimator.  In  those  regions,  the  window  is  overlapping 
two  range-independent  regions  with  with  wavenumber  values  very  closely  spaced. 
Examining  the  range-independent  regions  of  the  figure,  and  recalling  that  the  range 
axis  represents  ranges  at  the  center  of  the  aperture,  it  is  possible  to  determine  the 
range  of  data  included  in  each  aperture.  Thus,  the  result  at  range  =  1100  m  includes 
data  that  spans  the  range  interval  600  m  <  R  <  1600  m,  and  results  at  range  =  2100 
m  includes  data  that  spans  1600  m  <  R  <  2600  m.  By  including  all  the  data  spanned 
by  the  aperture  as  the  window  slides,  the  range-independent  regions  identified  from 
figure  3-14  can  be  extended  by  500  meters  at  each  end  to  clearly  demarcate  the 
range-independent  regions  for  TC3  as  shown  in  figure  3-15.  Using  this  method,  it 
was  determined  that  the  sediment  of  TC3  was  composed  of  a  background  region 
with  an  intrusion  consisting  of  a  rectangular  block  located  in  range  between  1100 
m  <  R  <  2900  m.  This  result  was  confirmed  by  the  workshop  organizers  when  the 
input  models  were  revealed  for  each  of  the  test  cases.  Having  identified  the  different 
regions,  wavenumber  estimates  were  made  using  the  entire  aperture  of  data  within  the 


90 


TC3  k  (r)  (50  Hz)  1000  m  AR 


Figure  3-15:  Wavenumber  estimates  for  range-dependent  sediment  of  Inversion  Work¬ 
shop  TC3  with  range-independent  segments  identified. 

3  range-independent  segments  as  shown  in  figure  3-16.  The  individual  wavenumbers 
estimated  for  the  intrusion  are  all  less  than  the  corresponding  values  estimated  for  the 
background  sediment  of  regions  1  and  3.  This  indicates  that  the  intrusion  consists 
of  a  harder  material  relative  to  its  surroundings  and  has  a  correspondingly  higher 
sound  speed.  Using  the  perturbative  inverse  method  previously  described,  this  result 
was  confirmed.  Figure  3-17  shows  the  inverted  sound  speed  profile  (dashed)  along 
with  the  starting  background  profile  (circles)  and  true  answer  (solid)  for  each  of  the 
regions. 

The  inversion  was  only  performed  at  a  single  frequency,  and  does  not  capture  the 
detailed  nature  of  the  true  profile,  but  matches  well  in  an  average  sense.  Resolution 
kernels  [65]  for  each  inversion  confirm  this  and  are  shown  in  figure  3-18  for  several 
depths  indicated  by  Zq.  In  general,  they  indicate  resolution  is  better  at  the  shallower 
depths  as  indicated  by  the  sharper  peaks.  Higher  resolution  could  be  obtained  by 
performing  wavenumber  estimates  and  inversions  at  higher  frequencies.  This  was 
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TC3  Spectra  Regions  1  and  3  (50  Hz)  TC3  Spectrum  Regions  2  (50  Hz) 


Figure  3-16:  Spectral  estimates  for  full  apertures  of  range-independent  regions  iden¬ 
tified  for  Inversion  Workshop  TC3. 


Seg.  1  &  3  Seg.  2 


Sound  Speed  (m/s)  Sound  Speed  (m/s) 


Figure  3-17:  Compressional  wave  speed  in  the  sediments  for  Inversion  Workshop 
TC3.  Open  circles  are  background  profile  used  for  perturbative  inversion,  solid  line 
is  ground  truth,  and  dashed  line  is  inferred  profile. 
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Intrusion 


Background 


Figure  3-18:  Resolution  Kernels  for  inversion  results  of  Inversion  Workshop  TC3. 
Background  represents  regions  1  and  3  and  Intrusion  represents  region  2. 

not  done  due  to  time  constraints,  and  because  the  data  were  not  provided  on  a  dense 
enough  range  grid  to  estimate  wavenumbers  without  aliasing.  Comparison  of  complex 
pressure  calculated  for  the  inverted  profile  matched  well  the  experimental  data  at  50 
Hz  and  100  Hz  and  shown  in  figures  3-19  and  3-20.  A  conclusion  from  this  test  case 
is  that  high-resolution  estimators  are  necessary  in  order  to  resolve  discrete  events 
that  occur  in  the  sediment  and  impact  on  the  complex  pressure  field.  It  is  clear  from 
the  analysis  in  figures  3-1— 3-3  that  the  intrusion  in  TC3  could  not  be  resolved  using 
classical  techniques. 
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Magnitude  of  Pressure  Field  (50  Hz) 


Phase  of  Pressure  Field  (50  Hz) 


Figure  3-19:  Complex  pressure  field  modeled  using  inferred  sediment  sound  speed 
profile  compared  to  inversion  workshop  TC3  data  at  50  Hz. 
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3.5  Wavenumber  estimation  in  random  media  -  In¬ 


ternal  Waves 

In  addition  to  making  wavenumber  estimates  for  waveguides  with  range  varying 
bathymetry  and  sediments,  measurements  in  realistic  ocean  environments  include 
the  effects  of  sound  speed  variability  in  the  watercolumn  due  to  temperature  fluc¬ 
tuations.  Although,  small  sound  speed  perturbations  may  be  indistinguishable  from 
other  causes  of  observation  noise  in  the  data  -  such  as  small  ranging  errors  or  range- 
dependence  in  the  sediments  or  bathymetry  due  to  sand  ripples  -  where  estimator 
characteristics  have  already  been  described,  larger  variations  must  be  considered. 
This  is  particularly  true  in  the  case  of  internal  waves  (IW),  where  it  is  well  known 
that  energy  transfer  can  occur  amongst  the  propagating  modes  [61].  The  temporal 
measure  of  sound  speed  variability  passing  due  to  internal  waves  is  on  the  order  of 
several  minutes  compared  to  a  few  seconds  for  a  propagating  acoustic  field  to  pass 
through  a  given  point  [14].  By  propagating  over  several  kilometers  through  multiple 
internal  wave  events,  the  received  field  takes  on  features  of  a  random  process.  In  a 
recent  paper,  Siderius  et.  al.  [74]  look  at  the  effect  of  sound  speed  fluctuations  on 
geoacoustic  inversion  results  for  broadband  acoustic  data  measured  on  vertical  line 
arrays.  In  the  wavenumber  analysis  discussed  for  this  thesis  and  the  experimental 
work,  the  effects  of  internal  waves  are  assumed  to  be  minimized  due  to  the  averaging 
effect  that  occurs  when  transforming  the  data  over  a  long  aperture  in  range  [14]. 
Further,  in  the  previous  discussion,  modal  evolution  was  discussed  in  the  context  of 
adiabatic  mode  theory  where  no  energy  is  transferred  between  modes  and  the  number 
of  propagating  modes  remains  constant.  However,  in  the  presence  of  internal  waves 
the  adiabatic  mode  assumption  is  violated  and  the  effects  on  local  wavenumber  es¬ 
timates  should  be  understood.  In  particular  it  is  of  interest  to  determine  if  internal 
wave  fields  cause  a  shift  in  the  estimated  wavenumbers,  if  they  cause  certain  modes 
not  to  be  resolved,  or  if  additional  modes  are  excited  and  can  be  detected. 
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To  understand  the  effects  of  internal  waves  on  the  forward  propagation  problem, 
a  recent  session  at  a  meeting  of  the  Acoustical  Society  of  America  (ASA)  [81]  focused 
on  benchmarking  acoustic  propagation  models  for  range-dependent  media  including 
internal  waves.  Although  focused  on  the  forward  problem,  these  studies  provided 
valuable  data  for  beginning  to  understand  the  wavenumber  estimation  problem. 

In  this  section,  synthetic  data  from  the  benchmark  session  is  used  for  estimating 
wavenumber  content  for  an  environment  both  with  and  without  internal  waves.  A 
description  of  this  model  and  other  range-dependent  cases  can  be  found  at  [75].  The 
waveguide  environment  was  a  Pekeris  model  with  a  flat  bottom  at  200  m  depth. 
Density  in  the  watercolumn  was  1  g/cm^  with  no  attenuation.  The  bottom  was 
modeled  as  a  fluid  with  sound  speed  of  1700  m/s,  density  1.5  g/cm^,  and  attenuation 
0.1  dB/A.  The  background  sound  speed  profile  was  generally  downward  refracting 
with  a  slight  duct  occurring  above  26  m  depth  and  was  given  by  [76] 


c{z)  =  1515  -f  O.OI62:,  z  <26  m 

c(z)  =  Co[l  -I-  a(e“*’  +  b  —  1)],  2:  >=  26  m. 


(3.43) 


where  Co  =  1490  m/s,  a  =  0.25,  and  b  =  (z  -  200)/500  [76]. 

To  the  background  sound  speed  was  added  a  perturbation  as  a  function  of  depth 
and  range  to  represent  the  internal  wave  field.  The  sound  speed  perturbation,  dc, 
was  given  by, 

5 

dc(z,r)  =  4{z/B)e~^^^ '^cos{Kir),  (3-44) 

1=1 

where  Ki  =  27r[2000  —  300(i  —  1)]“^  m“^  and  B  =  25  m.  Using  this  model,  the 
maximum  sound  speed  perturbation  is  about  7.5  m/s.  The  range  dependent  sound 
speed  profile  is  shown  in  figure  3-21.  Using  the  above  models,  PECan  [11],  was  used 
to  generate  the  acoustic  field  for  both  the  background  and  perturbed  sound  speed 
environments  from  0-20  km.  The  data  are  shown  in  figure  3-22  for  the  full  aperture 
and  for  5  km  sub-apertures.  The  difference  in  structure  of  the  modal  interference 
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Figure  3-21:  Range-dependent  sound  velocity  field  for  IW  wavenumber  estimates  [76]. 
Background  profile  is  on  left. 


patterns  observed  in  the  transmission  loss  plots  is  most  evident  at  long  ranges.  These 
data  were  generated  by  Gordon  Ebbesen  [21]  at  DREA  Canada  on  a  5  meter  range 
grid.  3000  m  sliding  window  AR  wavenumber  estimates  were  made  for  both  cases 
with  a  step  of  500  m  and  a  model  oder  of  200.  Additionally,  the  exact  mode  values  for 
the  unperturbed  model  environment  were  determined  using  KRAKEN.  The  resulting 
wavenumber  results  with  the  KRAKEN  results  overlayed  are  shown  in  figures  3-23 
and  3-24. 

The  number  of  modes  determined  by  KRAKEN  is  greater  then  the  number  esti¬ 
mated  using  the  sliding  window  method  indicating  that  not  all  the  modes  are  excited. 
Wavenumber  estimates  that  agree  with  the  KRAKEN  results  are  shown  by  a  dashed 
line.  The  plots  show  that  for  the  unperturbed  case,  6  modes  are  estimated  that  can 
be  correlated  with  the  KRAKEN  result.  For  the  IW  case,  8  modes  are  identified 
that  correlate  with  the  KRAKEN  results.  The  plots  show  the  wavenumbers  to  be 
very  stable  with  range,  indicating  range-independent  boundary  conditions.  This  is 
particularly  evident  for  the  the  higher  order  modes.  There  is  evidence  of  some  in¬ 
stability  of  the  estimates  for  the  IW  case  seen  as  undulations  about  the  predicted 
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Background  k^(r)  (100  Hz)  3000  m  AR 
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Figure  3-23:  Sliding  window  wavenumber  estimates  with  KRAKEN  results  (dashed) 
plotted  for  unperturbed  case. 
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Figure  3-24:  Sliding  window  wavenumber  estimates  with  KRAKEN  results  (dashed) 
plotted  for  IW  case. 
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value.  However,  this  result  might  be  expected  as  energy  is  being  exchanged  between 
the  modes  with  range  and  appears  as  a  scintillation  in  the  wavenumber  estimates 
with  range.  However,  there  is  no  apparent  shift  observed  for  the  modes  observed  in 
the  unperturbed  case  and  the  IW  case.  There  are,  however,  additional  modes  that 
are  energized  due  to  the  mode  coupling  which  are  consistent  with  those  predicted  by 
KRAKEN.  A  conclusion  from  this  study  is  that  for  the  environment  provided,  the 
internal  waves  couple  energy  into  modes  that  might  not  otherwise  be  excited  for  a 
given  source/ receiver  geometry  and  give  additional  information  that  could  be  used  as 
input  data  for  inversion.  This  is  a  good  result  for  the  relatively  weak  internal  wave 
field  provided.  However,  future  studies  could  include  an  examination  of  the  types  of 
limits  on  IW  perturbations  that  cause  these  observations  to  break  down.  A  second 
caveat  to  these  observations  relies  on  the  ability  of  an  estimator  to  resolve  individual 
modes.  The  exact  nature  of  the  coupling  and  the  wavelength  of  the  internal  wave 
structure  will  impact  the  ability  to  resolve  modes.  In  this  particular  case,  it  is  clear 
that  the  first  two  modes  predicted  by  KRAKEN  were  not  present  in  the  data.  How¬ 
ever,  the  first  two  modes  were  very  close  together  in  wavenumber  space  and  would 
require  a  long  aperture  to  resolve  them.  These  issues  are  less  troublesome  at  lower 
frequencies  and  shallower  depths  where  the  wavenumbers  are  less  densely  spaced.  In 
these  cases,  shorter  apertures  could  be  used  to  observe  the  same  effect  as  the  3  km 
apertures  used  in  this  study. 
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Chapter  4 


Modal  Mapping  Experiment 
(MOMAX) 

4.1  Introduction 

The  Modal  Mapping  Experiments  were  designed  as  an  efficient  means  for  measuring 
the  full  spatial  variability  of  a  point  source  acoustic  field  in  a  complex  shallow  water 
environment  [30].  Previous  experimental  efforts  have,  in  general,  employed  synthetic 
aperture  horizontal  arrays  generated  along  radials  to  measure  pressure  as  a  function 
of  range  at  a  constant  bearing.  Typically,  this  was  achieved  by  towing  a  hydrophone, 
suspended  at  a  fixed  depth,  away  from  a  moored  source,  or  alternatively,  towing  a 
source  away  from  a  fixed  receiver.  For  ranges  greater  than  a  few  wavelengths,  the  hor¬ 
izontal  wavenumber  spectrum  of  the  field  could  then  be  extracted  by  employing  the 
asymptotic  form  of  the  Hankel  transform.  These  efforts  were  effective  in  demonstrat¬ 
ing  the  range-dependent  nature  of  the  spectrum  due  to  range  dependent  waveguide 
properties.  However,  it  is  desirable  to  examine  the  effect  of  waveguide  properties  on 
the  spectrum  along  radials  at  all  angles  around  the  source  as  given  from  a  complete 
spatial  measure  of  the  propagating  field.  The  MOMAX  experiments  were  designed  to 
make  measurements  of  the  propagating  acoustic  field  for  a  point  source  on  a  spatial 
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Figure  4-1:  Sketch  representation  of  MOMAX  experimental  layout 

grid  given  in  absolute  geospatial  coordinates  of  latitude  and  longitude.  By  design, 
these  measurements  are  related  to  the  horizontal  wavenumber  spectrum  of  the  field 
by  application  of  the  appropriate  two-dimensional  transform  relation  as  discussed  in 
chapter  2.  An  overview  of  the  experimental  configuration  is  given  in  the  next  section 
followed  by  descriptions  of  specific  experiments. 

4.2  Experiment  Overview 

In  the  MOMAX  experiments,  data  are  acquired  on  an  array  formed  by  individual 
hydrophones  suspended  from  several  freely  drifting  buoys  equipped  with  precision 
GPS  navigation.  Figure  4-1  illustrates  the  operational  scenario  and  assets.  As  many 
as  four  buoys  might  be  drifting  at  a  given  time  and  synthetic  apertures  can  be  formed 
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using  data  acquired  on  a  single  buoy  or  several  buoys. 

For  the  usual  experimental  run,  the  acoustic  field  was  generated  by  a  single  pro¬ 
jector  operating  in  a  cw  mode  comprised  of  several  discrete  tones  between  20  and  475 
Hz.  At  any  given  time,  no  more  then  four  frequencies  were  transmitted  by  a  projec¬ 
tor.  During  one  experiment,  two  projectors  were  used  simultaneously  transmitting 
four  frequencies  each.  Several  different  acoustic  projectors  were  used  in  the  MOMAX 
experiments,  including  an  NRL  J-15-3  source  (50-475  Hz),  a  Webb  Research  organ 
pipe  type  source  (200,300  Hz),  and  an  IHI  source  (20-125  Hz).  Each  source  was  set  to 
operate  at  a  nominal  level  of  173  dB  re  IfiPa  at  the  set  frequencies.  Experiments  were 
run  with  the  source(s)  either  moored  on  the  seafloor,  or  suspended  from  the  ship  at 
fixed  depth  -  nominally  30  meters.  With  the  source  suspended  from  the  ship,  several 
different  operational  configurations  were  used  including  mooring  the  ship  to  achieve 
a  stationary  source  position,  allowing  the  ship  to  drift  freely,  or  towing  the  source 
with  the  ship.  These  configurations  play  an  important  role  in  the  interpretation  of 
the  data  as  will  be  shown  in  the  next  chapter. 

The  measurement  buoys  used  in  MOMAX,  schematically  depicted  in  Figure  4-2, 
were  developed  at  WHOI  specifically  for  these  experiments.  Each  buoy  contained 
a  calibrated  hydrophone  suspended  at  a  nominal  depth  of  30  meters  by  a  system 
designed  to  isolate  the  effects  of  wave  motion  on  the  vertical  motion  of  the  phone. 
Measurements  of  hydrophone  depth  using  pressure  sensors  showed  actual  hydrophone 
depths  to  be  greater  than  30  meters,  but  generally  stable,  keeping  depth  to  within  a 
meter.  A  typical  depth  record  for  a  single  experiment  is  plotted  in  figure  4-3.  The 
data  record  includes  data  taken  starting  before  deployment  and  continuing  through 
recovery  of  the  buoys.  A  depth  reading  of  zero  is  indicated  when  the  buoy  is  on 
deck.  Throughout  the  course  of  the  experiment,  it  was  observed  that  during  periods 
of  adverse  weather  at  the  surface,  hydrophone  depth  variation  increased  as  the  limits 
of  the  damping  system  were  exceeded.  It  should  be  noted  that  between  the  second 
MOMAX  experiment  and  the  SWAT  experiment,  the  hydrophone  suspension  and 
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Transponder,  8  lbs  wet 
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Upward  pointed  transducer 
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3dBBW  15-900  Hz 


Figure  4-2:  MOMAX  drifter  buoy  components. 
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SWAT  EXP  5  -  Hydrophone  Depth  (Larry) 


Figure  4-3:  Hydrophone  depth  record  for  typical  MOMAX  experimental  run. 

damping  system  was  modified.  Comparing  depth  records  from  both  experiments,  the 
system  used  on  the  original  buoys  was  superior  to  that  of  the  SWAT  experiment. 
This  observation  is  important  because  depth  data  was  not  available  for  the  MOMAX 
97  experiment  that  used  the  original  hydrophone  suspension  system. 

Acoustic  data  measured  at  the  phones  were  broadcast  back  to  the  ship  in  real¬ 
time  over  a  VHF  radio  link.  The  analog  signals  were  sampled  at  3255.2083333  Hz 
and  stored  to  computer  hard  disk.  The  buoys  were  also  equipped  with  GPS  receivers 
and  900  MHz  modems  for  transmitting  geographic  coordinates  back  to  the  ship  in 
real-time.  Using  a  mapping  program  developed  as  part  of  this  work,  a  real-time 
display  of  all  buoy  and  ship  positions  was  used  to  monitor  buoy  locations.  Additional 
GPS  receivers  were  used  to  measure  the  position  of  the  suspended  acoustic  sources. 
GPS  data  were  recorded  and  subsequently  processed  to  achieve  absolute  positioning 
accuracy  to  within  1  meter.  Further,  by  establishing  a  local  differential  GPS  system 
between  the  source  ship  and  each  buoy,  relative  positions  between  the  source  and 
each  receiver  were  measured  to  sub-meter  accuracy  [19]  [20]. 

Application  of  Hankel  transform  based  methods  to  the  measured  data  requires 
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accurate  measurements  of  the  acoustic  phase;  consequently  a  precise  time  base  was 
required  for  both  the  source  and  receiver  data.  The  source  signal  was  generated  by  a 
function  generator  with  an  external  clock  controlled  by  a  TRAK  GPS  clock  accurate 
to  within  200  nanoseconds  of  true  GPS  time.  The  same  clock  was  used  to  drive  the 
data  acquisition  system.  Synchronization  of  the  time  base  for  the  source  and  receiver 
systems  allowed  very  precise  phase  measurements  to  be  made  of  the  acoustic  signal. 

4.2.1  Data  Reduction 

The  acoustic  signal  was  generated  and  transmitted  as  a  sum  of  sinusoids  for  the 
frequencies  of  interest.  However,  the  Hankel  transform  method  described  is  for  a 
single  frequency.  To  look  at  a  single  frequency,  the  data  received  on  the  hydrophone 
were  quadrature  demodulated  at  each  of  the  transmitted  frequencies  [5]  by  taking 
the  FFT  of  approximately  1.25  seconds  of  data  with  time  determined  by  the  center 
of  the  FFT  window.  Complex  pressure  data  at  each  frequency  were  taken  as  the 
real  and  imaginary  parts  of  the  resulting  transform  at  the  appropriate  frequency 
bins.  Care  was  taken  to  ensure  that  the  resulting  spatial  sampling  was  on  the  order 
of  a  point  every  meter.  The  total  signal  is  treated  by  sliding  a  window  of  a  given 
number  of  points  and  performing  the  FFT  at  each  step.  A  half-window  overlap  was 
used  for  this  procedure.  Phase  accumulation  was  accounted  for  by  keeping  track 
of  time  accumulated  before  application  of  the  transform  at  each  step  and-  wrapping 
accordingly. 

Positioning  data  were  recorded  independently  on  a  separate  system.  The  acoustic 
and  navigation  data  were  merged  by  finding  the  union  of  times  where  both  acoustic 
and  position  data  were  recorded.  Positioning  data  were  acquired  every  second,  a 
slightly  faster  rate  then  the  demodulated  acoustic  data.  Consequently,  the  acoustic 
data  was  interpolated  onto  a  time  grid  common  to  the  GPS  data.  Linear  interpolation 
was  used  by  considering  the  real  and  imaginary  parts  of  the  complex  signal  separately. 
The  result  of  this  step  was  to  provide  data  on  an  xy-grid  relative  to  the  source 
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positioned  at  the  origin.  Similarly,  the  data  could  be  interpolated  onto  an  absolute 
position  grid  by  using  the  GPS  data  given  in  coordinates  of  latitude  and  longitude. 

4.3  MOMAX  97  and  SWAT  01 

Since  1997,  three  different  MOMAX  type  experiments  have  been  performed.  Two 
of  the  experiments  took  place  off  the  coast  of  New  Jersey  and  will  be  discussed  in 
this  thesis.  The  third  experiment  took  place  in  the  Gulf  of  Mexico,  about  145  nm 
west/northwest  of  Key  West,  Florida. 

4.3.1  Environment 

The  first  modal  mapping  experiment  took  place  in  about  70  meters  of  water  off  the 
NJ  coast  in  March  1997.  In  October  2001,  another  experiment  was  conducted  in  the 
same  area  under  the  Shallow  Water  Acoustic  Technology  (SWAT)  Memorandum  of 
Understanding  (MOU)  between  the  US  and  Japan.  The  overall  SWAT  experiment 
was  a  joint  effort  between  institutions  from  both  the  US  and  Japan.  In  addition  to 
the  MOMAX  group  from  WHOI,  other  US  participants  were  the  Naval  Research  Lab¬ 
oratory  (NRL),  and  the  University  of  Miami  (UM).  The  Japanese  participants  were 
under  the  direction  of  the  Technical  Research  and  Development  Institute  (TRDI). 

The  general  operational  area  for  both  experiments  was  within  the  STRATAFORM 
swath  mapping  survey  area  shown  in  figure  4-4.  STRATAFORM  is  a  program  of 
the  Office  of  Naval  Research  (ONR)  to  collect  high-resolution  geophysical  data  in 
shallow-water  regions  for  determining  continental  shelf  morphology  and  sediment 
properties  [53].  Presently  there  is  a  concerted  ONR  effort  [56]  by  both  geologists 
and  underwater  acousticians  to  characterize  this  region  in  terms  of  both  its  geologic 
and  acoustic  properties.  The  MOMAX  effort  is  concentrated  on  examining  the  re¬ 
lationship  between  local  sub-bottom  sediment  properties  and  low-frequency  acoustic 
propagation. 
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Specifically,  MOMAX  97  was  concentrated  in  two  regions  within  the  STRATAFORM 
area  centered  at  Ti.lTW  Longitude,  39.06°A^  Latitude  and  72.91°VF  Longitude, 
39.24°A’  Latitude  corresponding  to  the  HUNTEC  site  and  SWARM  site,  respectively. 
The  HUNTEC  site  was  where  an  extensive  high-resolution  3-D  seismic  survey  was 
conducted  with  the  discovery  of  a  complex  internal  sediment  structure  underlying  a 
more  benign  seafloor  surface  [16].  SWARM  refers  to  the  “  Shallow  Water  Acoustic  in 
Random  Medium”  experiment  that  took  place  in  1995  [3].  The  SWAT  experiment  re¬ 
visited  the  HUNTEC  and  SWARM  sites  with  focal  points  near  73.15°VU,  39.00°A^  and 
72.85°LF,  S9.30°N,  respectively.  Each  of  these  locations  is  within  a  region  referred  to 
as  the  “outer  shelf  sediment  wedge”  [16].  This  region  extends  along  the  edge  of  the 
New  Jersey  shelf  for  150  km  southwest  of  Hudson  Canyon.  The  entire  region  can  be 
characterized  by  a  layer  of  soft  sediments  of  varying  degrees  of  thickness  and  strat¬ 
ification  overlying  a  prominent  sub-surface  seismic  reflector,  designated  “R”.  In  the 
region  of  the  MOMAX  experiments  the  sediment  layer  is  thin  and  is  thought  to  be  a 
homogeneous  mixture  of  silty  and  sandy  clay  [16].  Additionally,  3-D  seismic  reflection 
surveys  have  identified  regions  within  the  wedge  where  the  sediment  layer  is  patchy. 

In  some  regions,  the  reflector  R  breaks  the  surface  and  sub-R  sediments  are  exposed. 

In  these  regions  a  complex  system  of  sub-marine  channels  has  been  identified  that  are 
believed  to  have  been  caused  by  an  erosional  episode  which  occurred  during  the  last 
(Wisconsinan)  period  of  low  sea  level.  Bathymetry  for  the  area  is  relatively  benign 
with  depth  gradually  increasing  along  lines  perpendicular  to  the  shelf  break. 

4.4  Summary  of  Measurements 
4.4.1  MOMAX  97 

During  the  course  of  MOMAX  97,  a  total  of  8  separate  experiments  were  conducted. 
Experiments  1-7  were  done  in  the  vicinity  of  HUNTEC,  and  experiment  8  was  con¬ 
ducted  at  the  SWARM  site.  Table  4.1  provides  a  summary  of  the  different  experiments 
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Exp.  No. 

Site 

Source 

Moe 

Curley 

Shemp 

1 

HUNTEC 

J-15 

no 

no 

YES 

2 

HUNTEC 

J-15 

no 

no 

YES 

3 

HUNTEC 

J-15 

YES 

YES 

no 

4 

HUNTEC 

J-15 

YES 

YES 

YES 

5 

HUNTEC 

WEBB 

YES 

YES 

YES 

6  ■ 

HUNTEC 

J-15 

YES 

YES 

YES 

7(a  &  b) 

HUNTEC 

WEBB 

YES 

YES 

YES 

8 

SWARM 

J-15 

YES 

YES 

YES 

Table  4.1:  Overview  of  source  and  buoy  deployments  for  MOMAX  97. 

including  site  location  along  with  the  type  of  source  and  the  buoys  deployed.  Source 
frequencies  for  the  J-15  were  set  at  50,  75,  125,  and  175  Hz  for  all  experiments,  and 
the  Webb  source  operated  at  200  and  300  Hz. 

Throughout  the  experiment,  environmental  data  were  collected  in  the  form  of  CTD 
data  for  determining  the  sound  velocity  in  the  water  column,  as  well  as  current  data. 
Each  buoy  and  the  source  were  also  equipped  with  temperature  sensors  that  recorded 
temperature  continuously  throughout  the  experiments.  As  shown  in  Figure  4-5,  sound 
velocity  in  the  water  column  was  generally  isovelocity  to  a  depth  of  between  30-35 
meters,  upward  refracting  between  35-60  meters,  and  isovelocity  for  depths  greater 
then  about  60  meters.  In  addition  to  the  water  column  measurements,  chirp  sonar 
data  were  collected  along  the  drift  tracks  of  the  individual  buoys  and  the  source  for 
many  of  the  experiments.  The  chirp  was  operated  for  a  frequency  range  of  3-6  kHz 
and  was  to  provide  detailed  bathymetry  data  over  the  course  of  each  experiment. 
Additionally,  the  chirp  data  could  be  used  to  indicate  qualitatively  the  stratigraphy 
in  the  region  of  the  individual  experiments.  Unfortunately,  for  this  experiment,  data 
recorded  on  DAT  tapes  for  archiving  were  not  usable  for  further  analysis  so  no  quan¬ 
titative  analysis  of  the  CHIRP  data  was  possible.  However,  as  part  of  the  SWAT 
experiment,  a  separate  CHIRP  survey  was  conducted  in  the  STRATAFORM  area  in 
the  vicinity  of  HUNTEC  and  SWARM.  As  this  data  becomes  available,  regions  of 
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Exp.  No. 

Site 

Source 

Source  Freqs. 

Moe 

Larry 

Shemp 

1 

SWARM 

IHI 

20,50,75,125, 

YES 

no 

YES 

J-15 

175,225,425,475 

2 

SWARM 

J-15 

50,  125,  475 

no 

no 

YES 

3 

HUNTEC 

J-15 

50,  75,  125,  175 

YES 

YES 

YES 

4 

HUNTEC 

J-15 

50,  75,  125 

YES 

YES 

YES 

5 

HUNTEC 

J-15 

50,  75,  125 

YES 

YES 

YES 

Table  4.2:  Overview  of  source  and  buoy  deployments  for  SWAT. 

overlap  with  the  MOMAX  97  experimental  data  will  be  sought  out  to  supplement 
the  environmental  characterization. 

Low-frequency  acoustic  data  was  taken  at  20,  50,  75,  and  125  Hz  with  a  source  level 
of  174  dB  re  1  //  Pa  using  the  J-15  source  suspended  from  the  ship.  It  was  intended 
that  the  ship  be  moored  or  at  anchor  to  provide  for  stationary  source  measurements. 
However,  due  to  weather  conditions  the  ship  could  not  be  anchored  and  the  exper¬ 
iments  were  conducted  with  the  ship  either  drifting  freely  or  towing  the  source  at 
a  constant  speed.  During  experiments  5  and  7  the  Webb  source  was  moored  and 
suspended  about  I  meter  off  the  seafloor  and  operated  at  200  and  300  Hz.  Unfortu¬ 
nately,  the  internal  clock  on  the  Web  source  became  unstable  during  the  experiments 
and  the  data  were  unusable. 


4.4.2  SWAT 

The  SWAT  experiment  was  comprised  of  5  different  experiment  runs,  two  near  the 
SWARM  site  and  3  near  HUNTEC.  The  source  and  buoy  deployment  overview  is 
given  in  table  4.2.  Frequencies  are  listed  separately  as  the  number  and  value  of  trans¬ 
mitting  frequencies  changed  for  each  run.  As  originally  planned,  each  of  the  SWAT 
experiments  was  to  use  two  sources  transmitted  4  frequencies  simultaneously.  Of 
particular  interest,  was  the  20  Hz  signal  provided  by  the  IHI  source.  Unfortunately, 
due  to  technical  problems,  the  IHI  source  was  only  operable  during  the  first  run.  As 
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Temporal  Sound  Speed  Variability  -  SWAT  EXP1 


-1  0123456789 

Time  (hrs)  m/s 


Figure  4-6:  Sound  speed  variability  during  SWAT  experiment  1  determined  from 
temperature  variability  recorded  on  NRL  T-string. 

in  MOMAX  97,  the  acoustic  experiments  were  supported  by  a  variety  of  environmen¬ 
tal  measurements  including  CTDs,  XCTDs,  and  XBTs,  along  with  temperature  and 
depth  sensors  on  both  sources  and  receivers.  Additional  assets  included  stationary 
temperature  string  (T-string)  moorings,  deployed  by  NRL,  that  spanned  the  water 
column  to  provide  a  continuous  record  of  the  temperature  variability  during  the  ex¬ 
periment.  Converting  the  temperature  variations  recorded  on  the  T-string  to  sound 
speed  using  Wilson’s  equation  [83],  assuming  a  constant  salinity  value,  an  example  of 
the  sound  speed  variability  measured  during  SWAT  experiment  1  is  shown  in  figure 
4-6.  Throughout  the  course  of  the  experiment,  the  temperature  profile  showed  the 
water  column  to  be  warm  near  the  surface  and  cooler  near  the  bottom  indicating  a 
downward  refracting  sound  speed  profile.  As  shown  in  the  picture,  the  location  of  the 
boundary  between  the  warm  and  cool  mixed  layers  was  time  dependent  with  a  vari- 
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ation  of  10  to  15  meters  in  depth  over  the  course  of  the  experiment.  This  variability 
would  be  of  concern  if  performing  wavenumber  transforms  over  the  entire  data  set. 
However,  it  will  be  shown  for  the  data  sets  considered  that  the  sound  speed  variability 
can  be  considered  locally  stationary. 

High-resolution  chirp  sonar  data  and  sediment  cores  were  also  performed  in  the 
area.  These  data  are  being  analyzed  separately  and  results  are  anticipated  for  future 
comparison  with  inversion  results  obtained  as  part  of  this  thesis  work.  In  the  following 
chapter,  data  from  the  MOMAX  97  and  SWAT  experiments  will  be  discussed  and 
analyzed. 
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Chapter  5 


Data  Analysis 

5.1  Introduction 

In  this  chapter  some  of  the  analysis  techniques  discussed  in  the  previous  chapters  are 
applied  to  the  MOMAX  experimental  data. 

5.2  Spatial  measurements  of  complex-pressure  fields 

The  MOMAX  experiments  were  novel  in  that  they  combined  measurements  of  acous¬ 
tic  field  data  with  navigation  data  provided  by  GPS.  Previous  synthetic  aperture  ex¬ 
periments  typically  used  either  acoustic  or  radar  ranging  systems  to  measure  source/receiver 
separation  distances  [48] .  The  nominal  accuracy  provided  by  these  systems  was  about 
one  meter,  where  radar  ranging  was  typically  better  at  longer  separation  distances 
and  acoustic  methods  used  for  near  ranges.  Using  GPS,  sub-meter  positioning  accu¬ 
racy  could  be  obtained  for  source/receiver  separation  distances  [19]  [20]  [31]  with  no 
restriction  on  the  ranges.  Further,  absolute  positioning  on  a  geospatial  grid  could  be 
determined  to  within  1  m. 

In  order  to  transform  the  complex  pressure  data  into  its  wavenumber  represen¬ 
tation,  it  is  necessary  to  accurately  measure  the  phase  of  the  signal.  Because  the 
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MOMAX  I  EXP2  SPATIAL  PRESSURE  FIELD  SHEMP  50  HZ 


-500 


-2500 


Meters  North  Meters  East 

Figure  5-1:  50  Hz  spatial  pressure  field  measured  with  source  located  at  the  origin 
of  the  coordinate  system.  Top  plot  is  pressure  magnitude,  center  plot  is  phase,  and 
bottom  plot  is  drifter  track  geometry. 


MOMAX  approach  is  new,  it  is  worthwhile  to  examine  the  phase  of  the  acoustic 
signal  in  order  to  both  comment  on  its  quality  and  also  possibly  gain  insight  into 
the  shallow  water  propagation  problem.  Figure  5-1  shows  50  Hz  complex  pressure 
measured  on  a  two-dimensional  grid  relative  to  the  source  position.  The  efiPects  of 
the  source/receiver  motion  are  clear  in  both  the  magnitude  and  phase  plots  as  the 
receiver  moves  out  in  range  from  the  source  and  then  turns  and  moves  broadside 
to  the  source.  Referring  back  to  the  modal  expressions  for  the  pressure  field  given 
by  equation  (2.36)  in  chapter  2,  the  phase  of  the  complex  pressure  field  measured 
in  space  is  proportional  to  range.  Using  the  radiation  condition  and  the  paraxial 
approximation  for  horizontal  wavenumbers,  given  by 


kn  =  ky/l  -  e„  k{l  -  e„/2)  for  Cn  <  1,  (5.1) 


117 


an  approximate  expression  can  be  derived  relating  the  time  rate  of  change  of  the 
phase  to  the  separation  range  rate  between  source  and  receiver  [33], 

dr  1  dcp 

_  _ 

dt  k  dt' 

where  A:  is  a  typical  wavenumber.  The  resulting  expression  can  be  integrated  in  time 
to  yield  a  relation  between  the  phase  and  the  source-receiver  range, 

r  =  To  +  p  (5.3) 

where  Vg  is  the  range  from  source  to  receiver  at  initial  time  tg.  From  these  relationships 
it  is  clear  that  as  source/receiver  separation  distances  increase,  phase  wraps  with  a 
positive  slope,  and  wraps  with  a  negative  slope  when  distances  are  decreasing.  These 
effects  are  evident  in  figure  5-2  where  complex  pressure  is  plotted  as  a  function  of 
time.  The  top  two  plots  are  complex  pressure  magnitude  and  phase,  and  the  bottom 
plots  are  separation  distance  between  source  and  receiver  measured  using  GPS  and 
calculated  using  equation  (5.3)  with  initial  range,  Vg,  assumed  known  as  determined 
by  a  point  GPS  measurement  at  tg.  Of  particular  interest  for  this  data  set  is  how 
well  the  phase  tracks  the  relative  source/receiver  positions  through  the  turning  point 
where  the  phase  reverses  itself  when  going  from  opening  in  range  to  closing  in  range. 
Figure  5-3  is  a  plot  of  range  rates  determined  from  both  GPS  measurements  and 
phase  measurements  at  50  Hz  using  equation  (5.2).  From  the  figures,  it  is  clear  that 
there  is  good  agreement  between  the  measurements  for  source/receiver  positioning 
using  both  GPS  and  acoustic  phase.  Using  this  type  of  analysis  on  data  acquired 
using  GPS  navigated  buoys  thus  shows  great  promise  for  doing  source  localization 
and  tracking  using  measured  phase  data  [58]. 
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GMT  Hours 

Figure  5-2:  50  Hz  data  plotted  as  a  function  of  time  with  corresponding  separation 
distance  between  source  and  receiver.  Initial  range  for  bottom  plot  determined  from 
a  single  GPS  measurement. 

5.3  Wavenumber  analysis  of  MOM  AX  data 

5.3.1  Geometrical  interpretation 

Given  the  previous  discussion  of  the  complex  pressure  field  measurements  from  the 
MOMAX  experiments,  it  is  of  interest  to  look  at  results  in  the  wavenumber  domain. 
As  discussed  in  chapter  2,  Hankel  transform  based  methods  are  to  be  used  for  examin¬ 
ing  horizontal  wavenumber  content  of  the  propagating  fields.  Because  these  methods 
are  based  on  a  radial  geometry,  the  measured  two-dimensional  spatial  pressure  data 
must  be  mapped  to  a  radial  grid.  Of  course,  this  is  not  ideal;  it  would  be  preferable 
to  use  the  full  two-dimensional  inverse  transform  relationship  given  by  (2.8),  but  the 
sparse  nature  of  the  measurements  precludes  its  use.  A  further  requirement  of  the 
Hankel  transform  based  methods  is  that  data  analysis  is  limited  to  data  where  the 
source  and  receiver  are  either  monotonically  opening  or  closing  in  range.  Given  this 
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EXP2  Shemp  Speed  vs  Time  (60  second  average  50  Hz) 


Figure  5-3;  Source/receiver  separation  speeds  determined  from  GPS  and  phase  mea¬ 
surements  of  50  Hz  complex-pressure  field. 

requirement,  spatial  data  must  first  be  mapped  as  a  function  of  range  indicating  the 
separation  distance  between  source  and  receiver.  From  this  mapping,  data  can  be 
selected  with  sufficient  range  apertures  for  resolving  individual  modes.  Referring  to 
figure  5-2,  data  measured  for  GMT  times  between  20  and  22.6  hours  is  monotoni- 
cally  increasing  with  range  and  suitable  for  Hankel  transform  based  processing.  The 
method  for  mapping  to  a  simple  range  grid  from  an  arbitrary  drift  path  is  illustrated 
conceptually  in  figure  5-4.  Having  mapped  the  synthetic  aperture  data  to  a  function 
of  range  relative  to  the  source,  other  considerations  must  be  given  to  interpretation 
of  results  given  the  full  two-dimensional  geometry  of  the  source/receiver  paths.  Much 
like  in  the  synthetic  aperture  experiments  for  linear  arrays,  assumptions  are  made  as 
to  the  local  nature  of  the  measurements.  Wavenumber  content  extracted  for  a  given 
range  interval  is  attributed  to  local  waveguide  properties  for  that  range  interval.  Sim¬ 
ilarly,  in  this  work,  local  waveguide  properties  are  attributed  to  measurements  made 
over  the  local  2-D  spatial  grid.  Thus,  local  inversion  results  for  a  given  sub-aperture 
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Figure  5-4;  Schematic  representing  mapping  of  spatial  pressure  data  to  a  range  grid. 

of  data,  must  be  given  in  terms  of  both  range  and  bearing  relative  to  the  source.  Al¬ 
though  not  much  more  will  be  said  regarding  this  issue  for  the  results  that  follow,  it 
becomes  an  issue  when  synthesizing  results  from  multiple  experiments  in  a  particular 
region.  Additionally,  when  mapping  the  data  to  the  range  grid,  consideration  must 
be  given  to  the  projection  angle  of  the  data  onto  the  local  radial.  Depending  on  the 
nature  of  the  spatial  variability  of  the  waveguide,  data  which  crosses  many  lines  of 
azimuth,  while  having  a  monotonic  range  separation,  may  have  to  be  interpreted  dif¬ 
ferently  then  data  which  is  taken  along  a  constant  bearing.  An  example  would  be  the 
case  where  the  bathymetry  approximates  a  wedge.  In  this  case  the  interpretation  of 
wavenumber  evolution  depends  on  the  orientation  of  the  receiver  track  to  the  source. 
For  a  track  where  the  receiver  travels  perpendicular  to  lines  of  constant  bathymetry, 
wavenumber  evolution  can  be  attributed  solely  to  range-dependence  of  the  environ¬ 
ment.  If  the  receiver  track  slowly  spirals  around  the  source,  wavenumber  evolution 
would  be  due  to  both  azimuthal-  and  range-dependent  changes  in  the  environment. 
For  this  thesis  work,  data  are  considered  for  cases  where  the  source/receiver  geometry 
is  dominated  by  their  separation  in  range  with  limited  changes  in  relative  bearing. 
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As  such,  when  mapping  the  data  to  a  range-grid,  azimuthal  symmetry  is  assumed 
over  the  apertures  considered.  In  future  work,  as  more  of  the  data  are  analyzed,  in¬ 
cluding  more  complicated  drift  paths,  more  sophisticated  techniques  may  have  to  be 
employed  to  exploit  the  full  spatial  measure  of  the  complex  pressure  field.  One  area 
to  explore  would  be  the  projection  theorems  commonly  used  in  image  processing.  As 
discussed  by  Lewitt  [46],  these  theorems  relate  a  one-dimensional  transform  along  a 
path  to  the  full  two-dimensional  transform  of  a  region  that  contains  the  path.  In¬ 
cluded  among  the  transform  theorems  is  the  projection-slice  theorem  [46]  that  relates 
a  one-dimensional  Fourier  transform  along  a  path  to  the  full  two-dimensional  Fourier 
transform.  This  may  be  a  way  to  exploit  the  sparse  nature  of  the  MOMAX  2-D 
pressure  field  measurements  to  extract  the  full  spatial  evolution  of  the  wavenumber 
spectra. 

5.3.2  Modal  Mapping  -  Single  Mode 

During  SWAT  01  measurements  were  made  using  a  20  Hz  source,  corresponding  to 
a  wavelength  of  about  75  m.  For  the  experimental  water  depths  of  approximately  70 
m,  and  a  sandy  bottom  with  sound  speed  1600  m/s,  theory  predicts  that  only  a  single 
mode  would  be  excited  for  a  Pekeris  waveguide  [28].  Figure  5-5  shows  the  source  and 
receiver  tracks  for  the  20  Hz  experiment,  where  the  source  location  is  indicated  by 
the  location  of  the  R/V  Endeavor. 
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SWAT/MOMAX  III  Experiment  1 


Figure  5-5:  Source/Receiver  geometry  for  SWAT  20  Hz  experiment. 


Data  are  considered  near  the  end  of  the  experiment  run  when  the  source  is  approxi¬ 
mately  stationary  and  the  buoys  are  between  7  and  9  km  from  the  source  and  drifting 
in  a  southerly  direction. 

The  single  mode  case  is  of  interest  because  the  modal  representation  of  the  pres¬ 
sure  field  is  given  by  a  single  term.  Using  the  adiabatic  approximation  for  a  range- 
dependent  waveguide,  modal  evolution  with  range  can  be  determined  without  trans¬ 
forming  to  the  wavenumber  domain.  The  adiabatic  mode  expression  for  the  pressure 
field  given  a  single  mode  is  written, 


P  = 


^^^exp  {i  ki{r')dr') 
\l ki{r')dr' 


(5.4) 


where. 


(5.5) 


Examining  the  phase  term,  <I>,  modal  evolution  of  the  single  mode,  ki{r),  can  be 
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determined  directly  by, 

h{r)  =  (5.6) 

This  method  was  applied  to  the  20  Hz  data  measured  on  the  buoys  Moe  and  Larry 
for  SWAT  experiment  1  as  described.  Pressure  magnitude  and  phase  measured  on 
Larry  are  shown  in  figure  5-6.  The  top  magnitude  plot  does  not  contain  the  deep  nulls 
typical  of  modal  interference  suggesting  the  interaction  of  several  modes.  Although 
the  magnitude  plot  appears  somewhat  noisy,  the  bottom  plot  shows  the  phase  to 
be  very  well  behaved.  Unwrapping  the  phase  data  and  applying  a  50  point  moving 
average  filter,  ki{r)  was  determined  using  equation  (5.6)  as  shown  in  figure  5-7.  Also 
shown  are  the  results  for  a  500  m  aperture  sliding  window  AR  transform.  Model 
order  selection  for  this  and  other  analysis  on  real  data  was  based  on  N/3  where  N  is 
the  number  of  data  samples.  The  number  of  points  for  discretization  in  wavenumber 
space  was  8192,  which  gives  a  bin  size  of  7.6699e-04  m~^,  which  is  greater  than  the 
variance  of  the  estimator  determined  by  the  method  of  Sakai  which  for  a  data  aperture 
of  500  m  and  an  assumed  SNR  of  10  dB,  is  approximately  1.25e-4  m“l.  However, 
even  given  the  low  resolution  in  wavenumber  space  for  the  AR  estimator,  along  with 
low  resolution  in  range  due  to  the  finite  step  side  of  the  estimator,  good  agreement  is 
shown  between  the  two  methods.  The  analysis  was  repeated  for  the  data  measured 
on  Moe  with  similar  results  as  shown  in  figures  5-8  and  5-9.  The  sliding  window 
transform  results  for  MOE  are  not  as  good  at  the  middle  ranges  because  there  was  a 
dropout  in  the  data  around  7900  meters  that  affects  the  AR  wavenumber  estimates 
as  the  transform  window  overlaps  the  region  of  missing  data. 
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Figure  5-6:  20  Hz  Complex  pressure  measured  on  Larry. 
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Figure  5-7:  Modal  evolution  for  single  mode  at  20  Hz  on  Larry. 
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Figure  5-8:  20  Hz  Complex  pressure  measured  on  Moe. 
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Figure  5-9:  Modal  evolution  for  single  mode  at  20  Hz  on  Moe. 
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In  this  section  results  are  presented  where  the  input  data  were  eigenvalues  determined 
from  the  MOMAX  measured  pressure  fields.  Range  independent  results  are  presented 
for  two  different  experiments  with  different  source/receiver  geometries.  For  each  case, 
eigenvalues  are  determined  and  used  as  input  data  to  a  perturbative  inverse  algorithm 
to  determine  an  effective  sound  speed  profile  for  the  sediment.  Magnitude  plots  of 
the  measured  data  are  compared  with  model  runs  using  the  profile  resulting  from  the 
inversion.  The  experiments  considered  are  from  the  MOMAX  97  data  set  measured 
in  the  STRATAFORM  region  near  New  Jersey. 

Figure  5-10  shows  the  source/receiver  configuration  for  MOMAX  97  experiment 
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MOMAX  I  EXP4  SHEMP  75  Hz 


Figure  5-11:  75  Hz  pressure  field  magnitude  for  MOMAX  97  experiment  4  plotted 
over  the  drift  track.  Drift  track  is  plotted  relative  to  source  position. 

4.  The  bathymetry  for  the  area  indicated  an  average  depth  of  about  77  meters.  For 
this  experiment,  the  receiver  is  considered  stationary  and  the  source  is  closing  in 
range  roughly  along  a  radial.  For  this  geometry  it  is  straightforward  to  apply  the 
Hankel  transform  based  methods  to  extract  wavenumber  content.  The  magnitude  of 
the  pressure  field  at  75  Hz  for  this  experiment  is  shown  in  figure  5-11  in  3-D  coor¬ 
dinates.  The  complex  pressure  as  a  function  of  range  is  shown  in  figure  5-12  where 
phase  is  shown  to  wrap  with  range.  The  aperture  for  this  experiment  was  only  about 
1300  m  so  a  single  transform  was  taken  to  determine  horizontal  wavenumber  con¬ 
tent.  The  spectrum  resulting  from  a  classical  spectral  estimator  using  Matlab’s  psd 
function  [60]  indicated  5  dominant  propagating  modes.  The  resulting  spectrum  is 
shown  in  5-13  where  the  result  was  interpolated  onto  a  finer  wavenumber  spacing 
grid  for  determining  the  exact  peak  locations.  Using  the  eigenvalues  determined  from 
the  peak  locations  of  the  spectrum,  a  background  model  was  determined  for  initial¬ 
izing  a  perturbative  inverse  algorithm  to  determine  compressional  wave  speed  in  the 
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MOMAX  I  EXP4  SHEMP  75  Hz  IP(t)l 
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Figure  5-12:  75  Hz  complex  pressure  plotted  as  a  function  of  range  for  MOMAX  97 
experiment  4. 
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Figure  5-13:  Horizontal  Wavenumber  Spectrum  for  full  aperture  75  Hz  data. 
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MOMAX  I  EXP4  SHEMP  75  Hz  INVERSION 
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Figure  5-14:  Starting  profile  and  inferred  profile  for  MOMAX  97  experiment  4.  Den¬ 
sity  was  1.8  g/cm^  and  no  attenuation  was  used. 

sediment.  The  model  was  determined  using  layers  with  sound  speeds  corresponding 
to  the  phase  speed  of  each  propagating  mode.  Given  a  number  of  layers  and  sound 
speeds,  The  basement  depth  and  corresponding  sound  speed,  are  adjusted  by  trial 
and  error  until  the  number  of  eigenvalues  for  the  model  matches  the  number  mea¬ 
sured  in  the  data.  As  the  basement  depth  is  increased  or  decreased,  each  layer  is 
stretched  or  compressed  accordingly.  KRAKEN  was  used  to  determine  the  model 
eigenvalues  for  the  background  model.  Using  the  perturbation  algorithm  described 
in  chapter  2,  the  background  model  was  updated  iteratively  until  convergence  was 
achieved  between  the  eigenvalues  of  the  updated  model  and  the  measured  values. 
The  background  starting  model  is  plotted  along  with  the  inferred  sound  speed  profile 
after  30  iterations  in  figure  5-14.  The  density  in  the  sediment  was  was  1.8  g/cm^ 
with  no  attenuation  in  the  model.  The  inferred  sound  speed  values  between  1500 
and  1800  m/s  for  the  upper  50  meters  of  sediment  are  consistent  with  sound  speeds 
for  sand  and  silty-sand  [36]  as  expected  for  the  STRATAFORM  region.  Using  the 
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MOMAX  I  EXP4  SHEMP  75  Hz 


Figure  5-15:  Comparison  of  75  Hz  measured  and  synthetic  pressure  magnitude  data 
for  inversion  of  MOMAX  97  experiment  4  . 

inferred  sound  speed  model  for  the  sediment,  along  with  measured  sound  speed  values 
for  the  water  column,  range-independent  synthetic  pressure  data  were  generated  and 
compared  to  the  measured  data  as  shown  in  figure  5-15.  The  inferred  model  predicts 
the  observed  modal  interference  pattern  well  with  the  exception  of  the  region  around 
1500  meters.  This  suggests  that  a  range  dependent  inversion  might  be  more  suitable 
for  this  region. 
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MOMAX  I  EXP2  SHEMP  TRACK 


Longitude  (degrees) 

Figure  5-16:  MOMAX  97  experiment  2  geometry  and  bathymetry.  Buoy  track  is 
shown  in  upper  left  and  source  track  in  bottom  right.  Water  depth  is  indicated  by 
the  colorbar  on  the  right  in  meters. 

The  source/receiver  geometry  and  bathymetry  for  a  second  experiment  from  MO¬ 
MAX  97  is  shown  in  figure  5-16.  For  this  experiment  both  the  source  and  receiver 
were  drifting  independently  along  different  paths.  This  geometry  is  a  representation 
of  that  described  by  Schmidt  and  Kuperman  [73]  where  the  source  and  receiver  are 
allowed  to  move  at  different  angles  relative  to  one  another  in  the  same  horizontal 
plane.  For  the  case  where  the  angles  between  the  velocity  vectors  of  the  source  and 
receiver  remain  constant,  the  radial  vector  between  the  source  and  receiver  can  be 
considered  constant  and  the  data  mapped  to  a  radial  grid  for  application  of  Hankel 
transform  based  methods.  The  data  represented  on  a  grid  with  its  origin  located  at 
the  source  is  shown  in  figure  5-17.  These  data  are  easily  mapped  to  a  range  grid  and 
the  full  aperture  spectrum  estimate  is  shown  in  figure  5-18. 

The  spectrum  determined  using  the  AR  spectral  estimator  showed  contributions 
from  6  dominant  modes  whose  eigenvalues  where  used  as  input  data  to  determine  a 
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MOMAX  I  EXP2  SHEMP  75  Hz 


Figure  5-17:  75  Hz  pressure  field  magnitude  for  MOMAX  97  experiment  2  plotted 
over  the  drift  track.  Drift  track  is  plotted  relative  to  source  position. 
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Figure  5-18:  Horizontal  wavenumber  spectrum  estimate  for  MOMAX  97  experiment 
2  data  at  75  Hz.  AR  estimate  is  solid  line,  PSD  estimate  is  dashed. 
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background  starting  profile  for  perturbative  inversion.  An  estimate  of  the  spectrum 
for  a  the  classical  estimator  with  a  Hanning  window  applied  to  the  data  is  shown  for 
comparison.  Of  particular  significance  is  the  lack  of  any  peak  in  the  classical  estimate 
at  a  wavenumber  around  0.295  m~^.  Initial  inversion  work  using  wavenumber  values 
from  the  dominant  peaks  in  the  classical  spectrum  gave  non-physical  results.  It  also 
serves  to  illustrate  the  subjective  nature  of  picking  which  peaks  are  significant  in 
the  spectral  plots  for  determination  of  wavenumbers  for  construction  of  background 
model.  Using  the  AR  results,  the  background  sound  speed  profile  and  the  resulting 
inferred  profile  are  shown  in  figure  5-19.  The  inferred  profile  shows  more  resolution 
because  the  layers  in  the  starting  profile  were  split  in  half  and  the  smoothness  con¬ 
straint  applied  to  the  least  squares  problem  of  the  perturbation  algorithm.  This  step 
represents  some  fine  tuning  to  the  general  inverse  algorithm  and  was  performed  with 
guidance  from  Dr.  Subramaniam  Rajan  [68].  Using  the  inferred  sound  speed  profile, 
synthetic  data  were  generated  using  KRAKEN  and  compared  to  the  measured  data 
as  shown  in  figure  5-20. 

The  model/data  comparison  is  in  general  very  good  with  the  modal  interference 
pattern  matching  very  well.  However,  in  order  to  match  the  model  to  the  data,  the 
ranges  for  the  data  had  to  be  shifted  by  370  meters.  It  is  unclear  if  this  shift  was 
necessary  to  accommodate  the  geometry  of  the  source/receiver  motion  or  if  there 
were  errors  made  in  the  measurements.  One  difficulty  of  the  moving  source  and 
receiver  problem  is  in  locating  the  result  of  the  inversion  geographically.  For  the 
moving  source/receiver  problem  it  is  not  possible  to  identify  a  specific  region  of  the 
waveguide  from  which  the  local  wavenumber  measurement  results.  This  would  be 
particularly  difficult  in  an  environment  that  has  properties  that  vary  in  both  range 
and  azimuth.  From  the  analysis  of  Hawker  [38]  and  Schmidt  and  Kuperman  [73], 
it  is  also  clear  that  the  source/receiver  motion  will  impart  a  shift  in  the  horizontal 
wavenumbers  of  the  the  spectrum  relative  to  that  of  a  stationary  experiment.  The 
shifted  wavenumbers  are  what  is  measured  and  the  inversion  is  based  on  using  the 
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MOMAX  I  EXP2  C(z)  Recovery 


Sound  Speed  (m/s) 

Figure  5-19:  Starting  profile  (thin)  and  inferred  (thick)  sediment  sound  speed  profile 
for  MOMAX  97  experiment  2.  Density  was  1.8  g/cm^  with  no  attenuation  used. 


TL  MOMAX  1  Exp  2  Shemp  75  Hz 


Figure  5-20:  Comparison  of  measured  and  synthetic  data  for  inversion  of  MOMAX 
97  experiment  2. 
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shifted  eigenvalues.  For  the  low  velocities  associated  with  the  drifting  source  and 
receiver  in  this  experiment  the  expected  wavenumber  shift  should  be  negligible.  This 
is  confirmed  by  the  good  agreement  between  the  measured  and  modeled  pressure  field 
magnitude  shown  in  figure  5-20  that  indicates  a  good  match  between  the  eigenvalues  of 
the  inferred  geoacoustic  model  and  the  measured  eigenvalues.  However,  these  results 
also  suggest  that  it  is  only  necessary  to  match  the  eigenvalues  of  the  propagating 
modes,  and  not  necessarily  the  propagation  environment,  to  generate  a  field  that 
matches  the  measured  data.  Thus,  when  determining  the  geoacoustic  model  for  a 
given  environment,  it  may  be  necessary  to  shift  the  measured  modal  eigenvalues  to 
account  for  source/receiver  geometry  and  motion  in  order  to  infer  the  true  geoacoustic 
properties  of  the  seabed. 
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SWAT  EXP2 


Longitude 

Figure  5-21:  Source  receiver  geometry  for  SWAT  experiment  2  with  region  of  source 
track  doubling  back  on  itself.  Segment  2  indicates  time  when  source  is  moving  away 
from  buoy  and  segment  3  indicated  time  when  source  is  closing  on  buoy. 

5.5  Modal  Mapping  -  Moving  Source  Analysis 

In  the  preceding  section,  a  discussion  was  given  for  the  moving  source/receiver  prob¬ 
lem  when  the  source  and  receiver  move  along  different  tracks  but  very  slowly.  In  this 
section,  the  problem  is  discussed  for  a  moving  source  that  traverses  the  same  space  in 
the  waveguide  twice.  In  chapter  2,  the  results  of  Hawker  [38]  and  also  Schmidt  and 
Kuperman  [73]  showed  that  a  shift  in  the  phase  of  the  complex  pressure  field  could  be 
expressed  as  a  shift  in  the  horizontal  wavenumbers  for  a  source  moving  at  a  constant 
velocity  where  the  shift  is  proportional  to  modal  group  velocity.  This  phenomenon 
was  observed  during  the  SWAT  experiment  when  wavenumber  estimates  made  for  a 
source  moving  out  and  back  over  the  same  waveguide  track  were  found  to  be  shifted 
relative  to  one  another.  The  experimental  geometry  for  this  experiment  is  shown  in 
figure  5-21  where  the  source  goes  out  along  a  radial  in  range  and  then  comes  back 
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SWAT  EXP2  (50  Hz) 


Meters  East 


Figure  5-22:  50  Hz  pressure  field  for  SWAT  experiment  2.  Source  is  moving  at  a 
constant  speed  of  2  m/s. 

again  along  roughly  the  same  track.  The  50  Hz  field  measured  for  this  run  is  shown 
in  figure  5-22.  The  pressure  field  magnitude  for  a  4  km  aperture  near  the  turnaround 
point  is  isolated  and  shown  in  figure  5-24.  Segment  2  refers  to  the  data  taken  along 
the  NW  leg  where  the  source/receiver  are  opening  in  range,  and  segment  3  refers  to 
the  return  leg  where  source/receiver  are  closing  in  range.  Sound  speed  variability 
measured  using  the  NRL  T-string  data  for  the  two  segments  are  shown  in  figure  5-23. 
For  both  segments,  the  location  of  the  mixed  layer  in  depth  and  the  temporal  variabil¬ 
ity  was  similar  as  determined  by  examining  the  mean  and  standard  deviations  of  the 
temperature  profiles.  The  modal  interference  pattern  for  the  two  tracks  are  very  sim¬ 
ilar  with  most  observable  differences  due  to  noise.  Using  a  measure  of  the  pressure  on 
the  hydrophone  for  a  period  of  time  when  the  50  Hz  source  was  turned  off,  SNR  was 
determined  to  be  about  20  dB.  Sliding  window  AR  transforms  with  1000  m  apertures 
were  applied  to  each  of  these  data  segments  for  estimating  horizontal  wavenumber 
content.  For  the  estimated  SNR  and  data  aperture  length,  the  theoretical  variance  of 
the  estimates  is  very  small  at  6.28e-7  giving  an  error  in  the  estimate  of  approxi- 
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Segment  2 


Segment  3 


Figure  5-23:  Sound  speed  variability  measured  on  NRL  T-strings  during  acoustic 
measurements  along  segments  2  and  3.  Colorbars  indicate  sound  velocity  in  m/s. 


mately  0.0003  m~^.  The  use  of  the  theoretical  variance  is  justified  by  the  observation 
that  the  AR  estimator  performed  as  well  or  better  then  predicted  on  the  time-series 
data  in  chapter  3.  Figures  5-25—5-28  show  the  wavenumber  evolution  for  the  two 
different  segments  using  two  diflferent  representations. 

Picking  off  the  peak  locations  for  the  different  mode  numbers  with  range  from 
the  waterfall  plots,  modal  evolution  with  range  are  shown  in  figure  5-29.  Difference 
between  maximum  and  minimum  wavenumber  estimates  determined  with  range  are 
an  order  of  magnitude  greater  than  the  estimated  measurement  error.  Thus,  even 
given  the  presence  of  noise  in  the  measured  signals,  the  wavenumber  estimates  are 
shifted  versions  of  one  another  for  the  overlapping  source  tracks.  This  is  consistent 
with  the  wavenumber  shift  for  a  moving  source  given  by  Hawker  [38]  and  also  ob¬ 
served  by  Rajan  et.  al.  [67].  From  this  observation,  a  method  for  correcting  for  the 
wavenumber  shift  is  suggested  along  with  a  direct  measure  of  modal  group  velocity 
from  the  data.  In  this  analysis,  the  broadband  effect  of  the  moving  source  is  removed 
by  analyzing  the  data  at  a  fixed  frequency.  The  resulting  wavenumber  shift  observed 
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Pressure  Magnitude  (50  Hz) 
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Figure  5-25:  Modal  evolution  for  50  Hz  data  of  segment  2. 
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SWAT  EXP2  (50  Hz)  Segment  3 


0.19  0.195  0.2  0.205  0.21  0.215 

k  (1/m) 

Figure  5-28:  Waterfall  plot  of  spectrum  amplitudes  with  range  for  50  Hz  data  of 
segment  3. 

in  the  data  due  to  source  motion,  as  predicted  by  both  Hawker  [38]  and  Schmidt 
and  Kuperman  [73]  through  the  phase  of  the  pressure  field,  must  be  accounted  for  in 
doing  geoacoustic  inversion. 

As  in  the  work  by  Rajan  et.  al.  [67],  the  starting  point  for  the  analysis  is  equation 
(2.40)  written  as 

K  =  ±  ~),  (5-7) 

'^ng 

where  v  is  the  source  speed,  positive  when  moving  toward  the  receiver,  and  the  source 
track  is  assumed  to  be  along  a  line  that  goes  through  the  receiver  location.  Estimates 
of  group  velocity  were  made  from  the  differences  obtained  from  the  observed  shift  in 
wavenumbers  determined  from  the  data.  For  the  individual  wavenumbers  determined 
along  each  track,  the  difference  in  wavenumbers  can  be  expressed  by, 

Afc;  =  *;(”*> =  (s.s) 

Vji  Vn 

where  v±  are  the  relative  velocities  for  the  different  tracks,  and  kn  is  the  eigenvalue 
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Modal  Evolution 


Figure  5-29:  Evolution  of  modes  1  and  2  with  range  showing  bias  due  to  source  motion 
over  segments  2  and  3. 


for  a  stationary  source.  For  a  source  moving  with  a  constant  velocity  back  and  forth 
along  the  same  track,  i.e.  =  — u_,  the  above  expression  can  be  simplified  to  give 
an  expression  for  the  group  velocity  as, 


2knV 


(5.9) 


In  the  above,  it  remains  to  estimate  kn  from  the  shifted  wavenumber.  For  a  source 
moving  out  and  back  along  the  same  path  at  a  constant  speed. 


k 


n 


(6.10) 


Combining  (5.8)  and  (5.10),  the  group  velocity  can  be  written  in  terms  of  sums  and 
differences  of  the  shifted  wavenumber. 


'^ng  — 


y*(v+)  u*{v-)  ^ 

rvn  f^n 


(5.11) 


Using  the  above  expression,  the  modal  group  velocity  can  be  calculated  from  the 
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shifted  wavenumber  observations.  The  wavenumber  estimates  corrected  using  range- 
averaged  group  velocities  for  the  apparent  Doppler  shift  are  shown  in  figure  5-30. 
The  expected  agreement  between  the  two  measurements  is  improved  and  an  average 
value  at  a  given  range  is  suggested  for  use  in  inversion  schemes.  The  calculated  modal 
group  velocities  of  784.4  m/s  and  668.9  m/s  for  the  first  two  modes  are  about  half 
the  expected  values  when  compared  with  the  water  sound  speed.  However,  the  source 
speed  of  2  m/s  was  very  slow  leading  to  small  observed  differences  in  the  wavenumber 
estimates  and  a  corresponding  lack  of  resolution.  Because  the  error  estimate  of  the 
expected  measure  was  small  compared  to  the  changes  in  individual  wavenumbers  with 
range,  some  range-dependence  in  the  medium  can  be  inferred.  Thus,  in  estimating 
the  error  in  the  measurement  of  the  modal  group  velocity,  the  error  in  the  sum  and 
difference  measurements  is  used,  where  the  expected  wavenumber  shift  with  range  is 
assumed  to  be  constant.  The  error  in  the  sum  of  the  wavenumbers  with  range  is  given 
by  and  taken  to  be  the  standard  deviation  of  the  measured  sum  as  a  function  of 
range  divided  by  two.  The  error  in  the  difference  is  given  by  and  taken  to  be  the 
standard  deviation  of  the  measured  wavenumber  differences  as  a  function  of  range. 
Using  the  above  definition,  the  error  in  the  group  velocity  measurement,  e„9,  can  be 
expressed  using  standard  techniques  as  a  percentage  of  the  expected  value  by  [79] 


Using  measured  values  in  the  above  expressions,  errors  in  the  group  velocity  measure¬ 
ment  due  to  wavenumber  estimation  errors  were  greater  than  40  %.  This  can  account 
for  much  of  the  error  in  the  observed  group  velocity  estimates.  Other  sources  of  error 
need  to  be  investigated,  including  the  effects  due  to  displacement  of  the  out  and  back 
tracks  for  this  data. 

Nevertheless,  motivated  by  these  experimental  results,  a  numerical  experiment 
was  conducted  using  modal  functions  and  eigenvalues  determined  using  KRAKEN 
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Modal  Evolution  (Corrected) 


Figure  5-30:  Doppler  corrected  wavenumber  estimates  for  modes  1  and  2  over  data 
segments  2  and  3. 


Depth  (m) 

Cp  (m/s) 

p  (g/cm^) 

0-75 

1510.0 

1.0 

75+ 

1600.0 

1.6 

Table  5.1:  Pekeris  waveguide  model  parameters  for  moving  source  experiment. 

for  a  Pekeris  waveguide  with  properties  given  in  table  5.1,  using  a  source  frequency 
of  50  Hz.  For  the  given  input  environment  and  source  frequency,  KRAKEN  also 
outputs  modal  group  velocities,  which  are  given  along  with  the  mode  numbers  and 
phase  velocities,  for  the  Pekeris  waveguide  in  table  5.2.  Using  the  mode  functions 
output  from  KRAKEN,  complex  pressure  fields  were  generated  as  a  function  of  range 
for  a  stationary  source,  and  for  sources  moving  with  velocities  of  ±  10  m/s.  The 
fields  were  generated  for  a  5  km  aperture  as  a  sum  over  the  modes  and  corresponding 
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Mode  No. 

kr  (1/m) 

Vp  (m/s) 

Vg  (m/s) 

1 

0.2055064180 

1528.71 

1499.07 

2 

0.1974793421 

1590.85 

1476.45 

Table  5.2:  Modal  eigenvalues  and  amplitudes  for  50  Hz  Pekeris. 


Mode  No. 

(1/m) 

(l/m) 

%  diff. 

„.true 

.9 

^jest 

.  .9  . 

%diff 

1 

0.2055064180 

0.2055054887 

4.5e-04 

1499.07 

1498.95 

0.008 

2 

0.1974793421 

0.1974895105 

5.1e-03 

1476.44 

1471.13 

0.365 

Table  5.3:  Wavenumber  and  group  velocities  estimated  from  moving  source  spectra. 


wavenumbers  using  the  simple  expression, 


p(r)  = 

n—\ 


^  ns^  nr 


K'T 


(5.13) 


where  2tnd  are  the  mode  amplitudes  at  the  source  and  receiver,  A:*  =  A:„(l  + 
with  Vng  the  modal  group  velocity,  v  is  the  source  speed,  and  w[r]  is  Gaussian 
white  noise.  The  complex  fields  with  SNR  of  40  dB  generated  for  the  different  source 
motions  are  shown  in  figure  5-31. 

Using  the  synthetic  data,  estimates  were  made  of  the  wavenumber  content  for  the 
different  source  velocities.  Figure  5-32  shows  the  spectral  estimates  for  the  10  m/s 
moving  source  fields.  The  two  estimates  are  shifted  from  one  another  with  respect 
to  the  stationary  value  of  the  wavenumbers  as  indicated  by  the  black  lines.  Using 
the  full  aperture  and  sampling  the  range  every  5  meters,  estimates  of  the  shifted 
wavenumbers  were  obtained.  Values  of  horizontal  wavenumbers  for  the  first  two 
modes  were  then  estimated  as  the  average  value  of  the  estimates  from  the  sources  of 
opposing  motion.  The  estimated  wavenumbers  compare  very  favorable  with  the  true 
values  for  the  given  environment  as  given  by  KRAKEN.  Group  velocities  were  also 
determined  which  also  compare  well  to  the  expected  values  as  shown  in  table  5.3.  The 
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Complex  Pressure  (50  Hz)  -  Pekeris 


Figure  5-31:  50  Hz  synthetic  pressure  fields  generated  for  moving  source  analysis. 
SNR  is  40  dB.  Source  speeds  were  10  m/s. 

results  from  this  experiment  suggest  that  the  wavenumber  differencing  technique  used 
to  estimate  discrete  wavenumbers  for  sources  moving  in  opposing  directions  over  the 
same  environment  is  quite  effective.  Estimates  of  wavenumbers  and  group  velocities 
agreed  well  with  the  true  values.  Relative  to  the  MOMAX  experiments,  that  had  a 
typical  source  velocity  of  2  m/s,  the  numerical  results  are  shown  for  a  high  source 
speed  in  order  to  emphasize  the  observed  wavenumber  shift  and  reduce  errors  in  the 
measurement  of  the  wavenumber  differences.  However,  the  numerical  experiment 
was  also  conducted  for  a  source  speed  of  2  m/s  and  yielded  a  small  degradation  in 
the  measurement  of  the  group  velocity  of  the  first  mode  which  was  estimated  as 
1478  m/s,  or  a  1.3  %  difference  from  the  true  value.  These  results  indicate  that 
the  wavenumber  differencing  technique  shows  promise  for  measuring  modal  group 
velocity  and  merits  further  study  both  numerically  and  experimentally.  In  particular. 
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Moving  Source  Spectra  (50  Hz) 


Figure  5-32:  Spectral  estimates  for  moving  source  data  with  stationary  values  of 
horizontal  wavenumbers  indicated  by  dashed  line.  Source  speed  was  10  m/s 

for  an  experimental  study,  careful  attention  should  be  given  to  reducing  errors  in 
the  wavenumber  measurements.  Much  could  be  done  in  this  regard  by  designing  an 
experiment  where  the  source  is  constrained  to  move  out  and  back  over  exactly  the 
same  range-independent  environment  where  the  water  column  is  well  mixed  and  the 
sound  speed  is  constant. 

5.6  Gelfand-Levitan  Results 

Before  leaving  the  analysis  chapter,  the  Gelfand-Levitan  method  is  considered  in  light 
of  some  of  the  practical  aspects  of  making  real  ocean  acoustic  measurements.  There 
are  three  issues  in  particular  that  affect  how  real  measurements  impact  the  applica¬ 
tion  of  the  Gelfand-Levitan  method.  The  first  issue  relates  to  the  assumption  of  no 
allowable  density  discontinuity  across  the  potential  layer  interfaces.  This  condition  is 
violated  in  almost  any  realistic  waveguide,  where  a  density  jump  typically  occurs  at 
the  water-bottom  interface.  It  will  be  shown  that  this  issue  can  be  resolved  through 
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application  of  an  impedance  matching  condition  applied  at  the  boundary.  The  second 
issue  is  related  to  the  integration  interval  in  wavenumber  space  of  the  point  source 
spectrum  to  determine  the  Fourier  transform  of  the  reflection  coefficient.  As  derived, 
the  integral  is  taken  over  all  wavenumber  space.  However,  in  an  experiment,  it  is  only 
possible  to  measure  real  propagation  angles  which  limits  the  wavenumber  integral  to 
±ko,  where  |fco|  is  the  minimum  water  wavenumber.  The  truncated  integration  inter¬ 
val  has  an  impact  on  the  transform  of  the  reffection  coefficient  similar  to  the  problems 
associated  with  short  aperture  spectral  analysis.  The  impact  of  truncating  the  inte¬ 
gration  interval  on  inversion  using  the  Gelfand-Levitan  approach  will  be  discussed.  A 
final  study  will  address  the  issue  of  determining  the  plane-wave  reflection  coefficient 
from  the  shallow-water  Green’s  function.  It  will  be  demonstrated  that  given  an  exact 
expression  for  the  Green’s  function,  the  reflection  coefficient  can  be  extracted  and 
sediment  properties  recovered. 

5.6.1  Effects  of  Density  and  Wavenumber  Interval 

A  numerical  study  was  conducted  to  examine  the  effect  of  a  density  discontinuity  at 
the  water-bottom  interface  on  inversion  using  the  Gelfand-Levitan  method.  A  deep 
water  problem  was  considered  with  a  sound  velocity  profile  in  the  sediment  consisting 
of  a  shallow  minimum  at  1  meter  in  depth,  followed  by  a  constant  positive  gradient 
with  depth.  The  profile  is  described  by, 

c{z)  =  1540  - (1540-  1515)^,  0<^<1 

c{z)  =  1515  +  0.97z,  1  <  z  <  145  (5.14) 

c(z)  =  1665,  z  >  145, 

and  is  shown  in  figure  5-35  as  the  exact  profile.  The  density  in  the  water  column 
was  assumed  to  be  1  gfcm^  and  in  the  sediment  to  be  1.5  g/cm^.  For  the  given 
environment,  the  reflection  coefficient  was  calculated  using  the  modified  Thomson- 
Haskell  propagator  matrix  approach  of  Mook  [52]. 
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Reflection  Coefficient 


Figure  5-33:  Reflection  Coefficient  at  20  Hz  with  and  without  density  discontinuity 
at  water-bottom  interface. 

To  address  the  density  discontinuity,  an  impedance  matching  condition  was  ap¬ 
plied  at  the  water-bottom  interface.  This  condition  is  presented  in  Tolstoy  and 
Clay  [80]  and  gives  a  means  for  expressing  the  reflection  coefficient  just  below  the 
interface  where  the  discontinuity  occurs. 

ry  _  Plfezo(l  T  R+)  Po^zl(f  ~  R+)  /g 

Pl^zo(l  +  R+)  +  Po^zl(l  “  R+) 

Here  is  the  reflection  coefficient  just  below  the  interface,  and  R+  is  the  reflection 
coefficient  just  above  the  interface.  The  reflection  coefficients  determined  using  the 
Thomson-Haskell  method  both  with  and  without  the  density  jump  are  shown  in 
figure  5-33.  The  effect  of  the  density  jump  is  evident  in  the  magnitude  plot  where  the 
reflection  coefficient  does  not  go  to  zero  for  large  k^.  Using  the  impedance  matching 
condition,  the  reflection  coefficient  determined  just  below  the  interface  matches  that 
for  no  discontinuity  and  approaches  zero  as  required  by  the  Gelfand-Levitan  method. 
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The  Fourier  transform  of  the  reflection  coefficient  as  a  real  function  of  depth 
and  is  plotted  in  figure  5-34  both  with  and  without  the  effects  of  the  density  jump. 
The  effect  of  the  jump  discontinuity  is  to  impart  an  extra  oscillation  in  the  result 
at  the  upper  depths  that  causes  an  offset  from  the  continuous  density  case  at  the 
deeper  depths.  The  Fourier  transform  results  were  used  as  input  data  to  the  Gelfand- 
Levitan  equation  to  invert  for  the  sound  velocity  in  the  sediment.  Figure  5-35  shows 
the  results  for  several  cases.  The  input  profile  is  shown  in  the  background  with  the 
Merab  [50]  result  plotted  directly  on  top  of  it.  This  result  is  from  the  corrected 
reflection  coefficient  shown  above  where  the  wavenumber  integral  to  determine  the 
Fourier  transform  was  taken  out  to  10  times  the  water  wavenumber.  The  highly 
oscillating  result  shows  the  effect  of  not  accounting  for  the  density  discontinuity.  The 
final  two  results  are  for  a  reflection  coeflScient  with  no  density  jump  and  with  the 
impedance  matching  condition  applied  where  the  wavenumber  integral  was  truncated 
at  ko-  This  result  shows  the  impact  of  truncating  the  wavenumber  integral  when 
taking  the  Fourier  transform  to  have  the  greatest  effect  at  large  depths. 

5.6.2  Shallow-water  Exact  Inversion 

The  final  numerical  experiment  was  simply  a  test  to  check  that  it  was  possible  de¬ 
termine  the  reflection  coefficient  from  the  shallow-water  Green’s  function  using  the 
relationship  given  by  (2.14).  It  was  thought  that  the  poles  in  the  Green’s  function 
might  make  the  expression  difficult  to  evaluate.  This  was  tested  using  the  shallow- 
water  Green’s  function  determined  for  the  Pekeris  waveguide  illustrated  in  figure 
5-36.  This  waveguide  was  chosen  from  an  example  given  in  Jensen  et  al.  [41]  with 
a  frequency  of  20  Hz,  source  depth  =  36  m,  receiver  depth  =  46  m,  and  sound  ve¬ 
locities  given  in  the  figure.  The  Green’s  function  was  determined  from  the  exact 
expression  for  the  reflection  coefficient  for  a  halfspace  and  compared  to  the  results 
of  OASES  [72]  and  KRAKEN  [59]  as  shown  in  figure  5-37.  The  agreement  was  good 
and  the  reflection  coefficient  was  determined  and  plotted  as  a  function  of  in  fig- 
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FFT{  R(k^) ) 


Figure  5-34:  Fourier  transform  of  reflection  coefficient  with  and  without  density  dis¬ 
continuity. 


Sediment  Sound  Speed  Reconstruction 


Figure  5-35;  Comparison  of  inferred  sound  speed  profile  in  sediment  using  the 
Gelfand-Levitan  method.  Results  are  shown  with  and  without  density  discontinu¬ 
ity  and  for  truncated  wavenumber  integrals. 
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Pekeris  Waveguide  Model 


Figure  5-36:  Pekeris  environment  used  to  calculate  shallow  water  Green’s  function. 

ure  5-38  where  the  density  jump  has  been  accounted  for.  The  Fourier  transform  of 
the  reflection  coefficient  and  the  resulting  inversion  are  shown  in  figures  5-39  and 
5-40.  Good  agreement  is  shown  between  the  inversion  and  the  input  model.  This 
demonstrates  that  a  measure  of  the  shallow-water  Green’s  function  might  be  used  for 
determining  the  reflection  coefficient  for  exact  inversion.  Work  is  in  process  for  using 
these  methods  on  Green’s  functions  determined  for  real  data  in  deep  water  [32].  The 
method  has  not  been  applied  to  shallow  water  data  at  this  point. 
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Horizontal  Wavenumber  (1/m) 

Figure  5-37:  Exact  Green’s  function  for  Pekeris  problem  compared  to  OASES  and 


KRAKEN. 
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Figure  5-38:  Reflection  coefficient  determined  from  shallow-water  Green’s  function. 


Half-Space  Waveguide  FFT  {  R  }  (20  Hz) 


Figure  5-39:  Fourier  transform  of  Pekeris  reflection  coeflBcient  with  effect  of  density 
discontinuity  removed. 


Half-Space  Waveguide  Inversion  Result  (20  Hz) 


Figure  5-40:  Sediment  sound  speed  proti?e  mf^rreS^ilsing  the  Gelfand-Levitan  method 
for  the  Pekeris  problem  in  figure  5-36. 
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Chapter  6 


Summary  and  Conclusions 


In  this  thesis,  the  relationship  between  the  spectrum  of  a  point  source  acoustic  field 
and  the  sediment  properties  of  a  shallow-water  waveguide  was  described.  Adiabatic 
mode  theory  was  then  used  as  a  guide  to  illustrate  how  modes  adapt  to  the  local  prop¬ 
erties  of  the  waveguide.  This  was  demonstrated  using  synthetic  data  for  a  waveguide 
with  properties  that  varied  in  three  dimensions.  The  resulting  modal  maps  could 
used  as  input  data  to  infer  local  waveguide  properties. 

As  part  of  the  ONR/SPAWAR  inversion  techniques  workshop,  synthetic  data  were 
provided  for  analysis.  For  TC3,  the  seafloor  was  flat  and  the  sediment  properties 
varied  with  range.  It  was  demonstrated  that  detection  of  range-variability  in  the 
sediment  using  wavenumber  analysis  required  high-resolution  methods  to  fully  resolve 
closely  spaced  wavenumbers  in  the  local  spectra.  A  high-resolution  method  based  on 
autoregressive  spectral  estimation  was  described  and  characterized  for  identifying 
local  modal  content  using  short-aperture  data.  The  high-resolution  estimator  was 
then  used  to  examine  wavenumber  content  for  a  waveguide  including  sound  speed 
perturbations  in  the  water  column  due  to  internal  waves. 

An  experimental  method  for  measuring  three-dimensional  acoustic  fields  was  de¬ 
scribed  for  use  in  extracting  spatially-varying  wavenumber  content.  The  resulting 
measurements  provided  data  on  synthetic  aperture  horizontal  arrays  of  arbitrary 
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shape  determined  by  the  source  and  receiver  motions.  The  full  2-D  horizontal  arrays 
were  sparsely  populated  requiring  that  data  be  mapped  to  a  range  grid  and  analyzed 
using  methods  based  on  the  zeroth-order  Hankel  transform.  Modal  eigenvalues  were 
extracted  from  select  pieces  of  data  and  used  as  input  data  for  perturbative  inversion. 
It  was  shown  that  interpretation  of  the  wavenumber  spectra  must  take  into  account 
source/ receiver  motions.  In  particular,  a  Doppler  shift  in  the  spectrum  was  observed 
for  an  experiment  where  the  source  was  moving  with  a  constant  velocity. 

In  addition  to  methods  and  analysis  based  on  perturbative  inversion  methods,  the 
Gelfand-Levitan  method  was  discussed  as  a  means  for  geoacoustic  inversion.  The  re¬ 
lationship  between  the  point-source  spectrum  and  the  plane- wave  reflection  coefficient 
of  the  water-bottom  sediment  was  discussed.  It  was  shown  that  the  full  spectrum, 
including  the  discrete  values,  could  be  used  to  determine  the  reflection  coefficient  for 
the  shallow-water  case.  A  simple  method  was  demonstrated  to  account  for  the  density 
jump  at  the  water-bottom  interface.  Finally,  the  effect  of  using  wavenumber  content 
corresponding  to  only  real  angles  of  propagation  was  examined.  This  work  suggests 
the  potential  use  of  MOMAX  type  data  for  doing  exact  geoacoustic  inversion. 

6.1  Conclusions 

Much  of  this  thesis  work  was  focused  on  a  means  to  provide  data  for  doing  pertur¬ 
bative  geoacoustic  inversion  for  spatially  varying  waveguides.  It  was  demonstrated 
that  a  tool  was  necessary  for  extracting  the  discrete  wavenumber  spectrum  for  short- 
aperture  data.  An  autoregressive  spectral  estimator  was  examined  for  use  as  a 
wavenumber  estimator.  It  was  demonstrated  that  the  estimator  performed  well  rel¬ 
ative  to  other  high-resolution  estimators  for  determining  the  frequency  content  of  a 
signal  with  noise.  The  same  analysis  was  applied  to  determine  wavenumber  content 
for  a  synthetic  acoustic  field  having  known  eigenvalues  with  similar  results.  An  ad¬ 
vantage  over  the  other  methods,  based  on  eigen-decomposition,  is  that  no  a  priori 
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information  regarding  the  signal  being  analyzed  is  required,  such  as  the  number  of 
propagating  modes.  For  the  range-dependent  inversion  workshop  TC3,  it  was  shown 
that  a  reduction  in  aperture  of  approximately  two  could  be  used  to  resolve  and  ex¬ 
tract  modal  content.  It  was  shown  that  the  location  and  shape  of  an  intrusion  in  the 
sediment  could  be  determined  exactly.  This  result  is  particularly  noteworthy  in  that 
for  the  inversion  workshop,  it  was  the  only  method  capable  of  determining  the  nature 
of  the  intrusion  in  the  sediment. 

In  the  study  of  internal  wave  effects  on  wavenumber  estimates,  it  was  observed 
that  the  sound  speed  perturbations  in  the  water  column  do  not  negatively  impact  the 
detection  of  propagating  modes.  Although  mode  coupling  is  observed  when  compared 
to  the  unperturbed  sound  speed,  the  effect  is  to  transfer  energy  to  modes  which  are 
consistent  with  the  boundary  conditions  of  the  waveguide.  The  net  effect  appears 
to  be  that  mode  coupling  yields  more  information  which  might  be  used  in  the  per¬ 
turbative  inverse  algorithm,  i.e.,  modes  that  may  not  otherwise  be  excited  are  now 
identified. 

For  the  MOMAX  experiments  described,  it  was  shown  that  the  synthetic  apertures 
created  by  the  drifting  GPS  navigated  buoys  were  of  very  high  quality.  Relative  source 
and  receiver  ranges  and  range  rates  could  be  determined  directly  from  the  phase 
measurements.  For  the  single-mode  experiment,  agreement  was  obtained  between  the 
AR  estimator  and  the  phase  derivative  for  examining  the  range  evolution  of  horizontal 
wavenumber  for  the  propagating  mode.  For  a  case  where  the  aperture  opened  along 
a  radial,  wavenumber  estimates  were  used  as  input  data  for  perturbative  inversion. 
The  resulting  sediment  sound  velocity  profile  was  consistent  with  a  sandy  sediment 
as  expected  for  the  region.  Wavenumber-based  inversion  was  also  performed  for  a 
moving  source/receiver  configuration.  The  field  determined  from  the  profile  resulting 
from  the  extracted  wavenumbers  matched  the  data  quite  well.  Although  the  moving 
source/receiver  configuration  made  it  difficult  to  assign  the  inversion  results  to  a 
specific  geographic  location,  the  results  suggest  that  wavenumber  analysis  techniques 
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as  a  means  for  generating  transmission  loss  curves  for  rapid  assessment  in  situations 
where  source  and  receiver  positions  are  in  constant  motion.  Related  to  this  result  was 
the  observation  that  wavenumber  estimates  where  shifted  for  an  experiment  where  the 
source  traversed  the  same  local  environment  along  an  outgoing  and  incoming  track. 
It  was  shown  that  a  better  estimate  of  the  true  wavenumbers  could  be  obtained  by 
correcting  for  the  Doppler  shift  in  the  spectrum  due  to  the  source  motion.  From  this 
analysis,  a  method  was  suggested  from  which  modal  group  velocity  could  be  estimated 
from  the  observed  Doppler  shift  for  a  moving  source.  The  results  of  a  numerical  study 
based  on  wavenumber  estimates  for  a  modal  field  with  shifted  eigenvalues  verified 
these  results  and  suggests  further  experimental  work. 

Finally,  the  numerical  studies  using  the  Gelfand-Levitan  method  suggest  potential 
use  of  MOMAX  type  data  for  performing  exact  geoacoustic  inversions.  It  was  shown 
that  for  wavenumbers  corresponding  to  real  angles  of  propagation,  inversion  results 
could  be  obtained  to  a  depth  of  one  or  two  wavelengths  in  the  sediment.  It  was 
shown  that  the  density  jump  at  the  water-bottom  interface  could  be  accommodated 
by  requiring  a  single  point  measurement  of  the  density  at  the  interface.  It  was  also 
shown  numerically  that  the  shallow- water,  depth-dependent  Green’s  function  could 
be  used  to  determine  the  plane-wave  refiection  coefficient  at  the  bottom  for  use  in 
the  Gelfand-Levitan  equation. 

6.2  Suggestions  for  Future  Work 

The  autoregressive  wavenumber  estimator  used  for  this  work  provides  reliable  wavenum¬ 
ber  estimates  only  for  regions  of  local  waveguide  invariability.  In  the  case  of  a  constant 
bathymetric  slope  where  the  sediment  properties  follow  the  bathymetry,  the  require¬ 
ment  of  local  invariability  is  a  major  deficiency.  One  way  wavenumber  estimates 
might  be  improved  in  this  case  is  through  the  implementation  of  a  time-varying  au¬ 
toregressive  estimator  (TVAR).  The  TVAR  estimator  allows  for  the  coefficients  of 
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the  parametric  signal  model  to  be  time- varying,  or  in  this  application,  range- varying. 
The  effects  of  internal  waves  on  wavenumber  estimates  is  another  area  that  deserves 
further  work.  The  limits  of  wavenumber  estimation  in  the  presence  of  stronger  and 
stronger  internal  wave  fields  should  be  studied. 

A  related  area  for  research  is  the  effect  on  wavenumber  estimates  due  to  the 
shape  of  intrusions  in  the  sediment  for  a  waveguide  with  constant  bathymetry.  For 
the  inversion  workshop,  a  rectangular-shaped  intrusion  was  identified.  An  interesting 
study  would  be  the  identification  of  other  regular  shaped  intrusions  such  as  a  semi¬ 
circle  or  triangle.  This  might  complement  work  relating  to  the  recent  GeoClutter  [56] 
experiments,  where  discrete  features  in  the  sediment  are  studied  as  contributors  to 
both  propagating  and  scattered  acoustic  fields.  This  work  might  also  provide  a  means 
to  determine  range  dependent  features  of  the  environment  which  could  be  used  along 
with  the  global  search  methods  for  doing  geoacoustic  inversion. 

From  the  experimental  work,  it  was  shown  that  the  measured  spectra  are  influ¬ 
enced  by  the  source  and  receiver  motions.  A  method  was  described  for  inferring 
the  group  velocity  from  the  shifted  wavenumber  spectra.  An  experiment  designed 
to  specifically  measure  shifted  spectra  for  group  velocity  measurements  is  suggested 
where  the  source  is  towed  at  a  high  rate  of  speed  to  maximize  the  Doppler  shift  and 
increase  the  observed  wavenumber  differences  for  increased  fidelity  in  the  difference 
measurements.  A  similar  experiment  might  also  be  designed  to  provide  data  along  a 
radial  at  constant  bearing  for  a  long  synthetic  aperture  to  provide  data  for  estimating 
the  plane-wave  reflection  coefficient  from  the  shallow-water  Green’s  function.  This 
would  require  additional  work  in  the  estimation  of  the  Green’s  function  from  data.  In 
this  thesis  work,  the  far-field  approximation  to  the  Hankel  function  was  used,  where 
an  exact  Hankel  transform  should  be  employed  for  the  best  estimate  of  the  total 
spectrum. 
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